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1.  INTRODUCTION 

Structures  designed  as  venting  enclosures  have  been  proposed 
for  use  in  the  containment  and  suppression  of  the  explosions  result- 
ing from  the  accidental  detonation  of  explosives  during  their 
manufacture,  handling,  and  storage.  These  structures  must  provide 
total  containment  of  fragments,  partial  suppression  of  the  shock 
wave,  and  a controlled  release  of  the  internal  quasi-static  pressure. 
Because  of  the  relatively  low  frequency  of  accidental  detonations, 
the  most  cost  effective  approach  to  the  design  of  such  structures 
appears  to  be  one  utilizing  inexpensive  and  easily  erectable 
structures  which  can  be  replaced  after  an  explosion. 

There  are  several  approaches  which  one  could  take  in  designing 
structures  to  contain  the  internal  explosion  while  undergoing  large 
plastic  deformation.  The  most  sophisticated  approach  would  involve 
the  use  of  a large  scale  finite  element  computer  code  to  model  the 
complete  structure  and  solve  for  the  dynamic  response  using  an 
elastic-plastic,  strain-rate  dependent  material  model.  In  addition 
to  being  beyond  the  capabilities  of  most  structural  analysis  computer 
codes,  such  analyses  are  very  expensive  and  time-consuming,  espe- 
cially when  being  used  in  a design  iteration  mode.  At  the  other 
end  of  the  spectrum,  a designer  could  utilize  simple  handbook  type 
of  equations  to  analyze  individual  components  of  the  structure 
independently.  Such  a design  procedure  would  be  very  quick  and 
inexpensive;  however,  the  resulting  design  would  be  highly  conserva- 
tive and  would  involve  a large  weight  and  cost  penalty  in  its 
construction.  The  objective  of  the  current  work  is  to  provide  the 
designer  with  easy-to-use  design  tools  which  will  enable  him  to 
simplify  the  design  procedure  without  paying  a large  weight  and  cost 
penalty. 

To  be  weight  and  cost  effective,  structures  should  be  designed 
so  as  to  undergo  large  plastic  deformations  under  the  loading  due  to 
an  accidental  explosion.  In  designing  structures  for  large  deflec- 
tion, dynamic,  plastic  response,  a useful  simplification  is  to 
neglect  elastic  effects  and  treat  the  material  as  a rigid-plastic 
material.  It  will  be  demonstrated  that  theories  based  on  such 
simplifications  can  provide  highly  accurate  design  procedures  for 
designing  suppressive  shielding  structures. 
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The  work  presented  in  this  report  is  separable  into  three  distinct 
tasks : 


(1)  Development  of  the  above  mentioned  structural  response 
model  that  is  applicable  to  suppressive  shield  design  (Chapter  2,  3, 
and  4;  Appendices  A and  B) . 

(2)  Development  of  a hydrocode  method  for  generating  the 
complete  flow  field  for  shock  interaction  with  suppressive  shield- 
type  panels.  These  calculations  provide  time-dependent  blast  loading 
transmission  functions  (Chapter  5)  for  use  in  structural  response 
calculations.  Load  functions  determined  with  this  method  are  used 
with  the  structural  response  model  in  a sample  calculation  for  a 
1/4-scale  category  1 suppressive  shield. 

(3)  Development  of  a model  for  computing  time-dependent 
internal  pressure  loads  produced  by  burning  propellant  and  pyrotech- 
nics inside  a vented  enclosure  (Chapter  6) . Calculations  for  a 
category  5 suppressive  shield  are  presented. 
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2.  BACKGROUND  FOR  STRUCTURAL  RESPONSE  MODEL 

In  formulating  approximate  theories  for  the  behavior  of  structural 
systems,  three  basic  criteria  must  be  addressed.  First,  an  equation 
which  expresses  an  equilibrium  balance  must  be  stated.  In  order  to 
apply  this  equilibrium  equation,  it  is  necessary  to  invoke  boundary 
conditions  which  constrain  the  possible  stress  and  displacement 
fields.  Problems  involving  dynamic,  plastic  deformation  of  common 
structural  elements  are  generally  much  too  complex  to  solve  analy- 
tically and  are  beyond  the  capabilities  of  most  structural  analysis 
computer  codes.  In  developing  approximate  solutions  for  such 
problems,  it  is  essential  to  have  a clear  formulation  in  which  the 
influences  of  the  approximations  can  be  clearly  assessed.  ; 

A dynamic  equilibrium  equation  which  is  very  useful  for 
formulating  approximate  solutions  in  plasticity  is  the  principle 
of  virtual  velocities  (reference  1) 

Js  pi  flids-  Jv  muiuidv  = Jv  Qj  <5jdv  (1) 

where  ui  and  <5j  are  kinematically  admissible  velocity  and  generalized 
strain  rate  fields.  The  generalized  stress  field  Qj  and  the  accelera- 
tion field  u i are  statically  admissible.  Statically  admissible  fields 
are  those  which  satisfy  the  equilibrium  constraints  of  forces  and 
accelerations,  while  kinematically  admissible  fields  are  those  which 
satisfy  the  kinematic  boundary  conditions  and  constraints  of  the 
problem.  More  comprehensive  discussion  of  these  features  can  be 
found  in  Martin  (reference  2)  and  Save  and  Massonnet  (reference  3) . 
Equation  (1)  is  a statement  of  instantaneous  power  balance.  The  first 
integral  represents  the  rate  at  which  the  external  forces  are  doing 
work.  The  power  of  the  inertial  forces  is  expressed  by  the  second 


1.  Martin,  J.  B. , "Mode  Approximations  for  Impulsively  Loaded 
Structures  in  the  Inelastic  Range,"  in  Structures,  Solid 
Mechanics,  and  Engineering  Design,  edited  by  M.  Te'eni,  Wiley 
Interscience , pp.  1227-1248  (1969) . 

2.  Martin,  J.  B. , Plasticity,  Fundamentals  and  General  Results,  The 
MIT  Press  (1975) . 

3.  Save,  M.  A.  and  Massonnet,  C.  E.,  Plastic  Analysis  and  Design  of 
Plates,  Shells  and  Disks,  Nor th-Holland  Publishing  Company  (1972) . 
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integral  or  alternately  this  can  be  interpreted  as  the  rate  of  change 
of  kinetic  energy.  The  last  integral  represents  the  instantaneous 
rate  at  which  energy  is  being  dissipated  in  plastic  work.  It  is  impor- 
tant to  emphasize  that  the  set  of  &i  and  is  independent  of  the 
set  of  Qj  and  ii^.  Formally,  we  may  say  that  the  kinematically 
admissible  field  is  independent  of  the  dynamically  admissible  field. 

If  the  two  fields  are  associated,  that  is  &i  are  the  actual  veloci- 
ties caused  by  the  loads  Pi,  then  Equation  (1)  reduces  to  a direct 
energy  balance.  However,  the  great  utility  of  the  principle  of 
virtual  velocities  is  that  it  does  hold  when  the  fields  are  not 
associated. 

A useful  result  for  assessing  the  errors  in  approximate  solutions 
caused  by  relaxing  the  admissibility  restrictions  on  the  displacement 
and  stress  fields  is  Hill's  principle  of  maximum  plastic  work 
(reference  4)  which  states  that  the  true  increment  in  plastic  work 
required  to  bring  about  a given  set  of  strain  increments  is  greater 
than  the  increment  in  plastic  work  which  would  be  required  to  bring 
about  this  same  strain  increment  by  any  other  distribution  of  stresses 
which  conforms  to  the  same  yield  criterion.  Hence,  for  any  kinemat- 
ically admissible  displacement  field,  in  order  to  dissipate  a given 
amount  of  energy  as  plastic  work,  an  approximate  stress  field  which 
does  not  satisfy  all  of  the  requirements  of  the  true  stress  field 
will  cause  larger  strains,  and  hence  larger  displacements  than  the 
true  stress  field.  While  this  result  offers  no  measure  as  to  the 
accuracy  of  any  approximation,  it  does  guarantee  that  the  approximate 
solution  will  overpredict  displacements,  a conservative  design 
approach. 

For  a rigid-plastic  material,  the  Levy-Mises  equations  provide 
the  general  relationship  between  the  components  of  stress  and  the 
plastic  strain  increments. 


dEij  = dX  °ij 


(2) 


If  the  von-Mises  yield  criterion. 


[(a  -a  )2  + (a  -a  ) 2 + (a  -a  )2 
L'  xx  yy'  ' yy  zz'  ' zz  xx' 


+ 6 (a  2+a  2+a  2)]  = 2a  2 

' xy  yz  zx  o 


(3) 


is  used,  then  the  constant  dX  is  found  to  be  proportional  to  the 
increment  in  effective  plastic  strain.  The  components  of  the 
stress  deviator  tensor  are. 


4.  Hill,  R. , The  Mathematical  Theory  of  Plasticity,  Oxford 
University  Press  (1956)  . 
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(o 


a 


11 


ij  am5 6ij 

(4) 

+a22+033)//3 

(5) 

Physically,  the  Levy-Mises  equations  require  that  the  increments 
of  strain  be  orthogonal  to  the  stress  deviator.  While  it  is  often 
useful  to  relax  this  requirement  slightly  in  order  to  obtain 
approximate  solutions,  its  influence  should  always  be  considered 
in  selecting  admissible  fields  of  generalized  stress  and  strain. 


An  excellent  review  of  the  range  of  validity  for  the  use  of 
rigid-plastic  theories  in  structural  dynamics  was  given  by  Symonds 
(reference  5) . One  of  the  most  successful  applications  of  limit 
analysis  to  rigid  plastic  structural  elements  has  been  the  use  of 
stationary  mode  solutions  (reference  6) . Martin  (reference  1)  has 
given  several  results  useful  for  determining  the  "best"  stationary 
mode.  However,  his  results  are  applicable  only  to  impulsively  loaded 
linear  problems  in  which  no  interaction  between  components  of  the 
generalized  stress  vector,  Q j , occurs.  Very  little  analytical  work 
has  been  done  on  problems  involving  finite  duration  loading  and 
interacting  dissipation  mechanisms. 

For  any  particular  problem,  once  a stationary  mode  of  the  form 


Ui(x,t)  = ♦■(x)W(t) 


(6) 


has  been  chosen,  a suitable  set  of  generalized  strain  rates 


<$j(x,t)  = \p.(x)¥(W(t))  (7) 

consistent  with  the  assumed  mode  can  be  determined.  The  generalized 
stresses  are  chosen  in  a manner  such  that  Q.d^  is  a rate  of  energy 
dissipation  and  we  assume  -* 


Q(x,t)  = ct  (x)  2 (W( t) ) . 


(8) 


5.  Symonds,  P.  S. , "Survey  of  Methods  of  Analysis  for  Plastic 

Deformation  of  Structures  Under  Dynamic  Loading,"  Technical  Report 
BU/NSRDC/1-67 , Division  of  Engineering,  Brown  University  (1967). 

* 6.  Martin,  J.  B. , and  Symonds,  P.  S.,  "Mode  Approximations  for 

Impulsively-Loaded  Rigid-Plastic  Structures,"  Journal  of  the 
Engineering  Mechanics  Division,  ASCE,  Vol.  92,  EM5,  43, 

* pp.  43-66  (1966). 
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W is  the  displacement  of  a reference  point  in  the  structural  element. 

The  form  chosen  for  the  generalized  stresses  implies  they  are  inde- 
pendent of  the  generalized  strain  rates.  For  many  materials  this  *■ 

is  a valid  assumption?  however,  most  common  grades  of  structural 
steel  show  a very  strong  influence  of  strain  rate  on  the  value  of 
the  uniaxial  yield  stress.  Bodner  (reference  7)  has  assessed  the  ? 

significance  of  strain  rate  effects  in  solving  dynamic  plastic 
problems.  He  concluded  that  for  structures  with  relatively  simple 
dissipation  modes,  such  as  a cantilever  beam,  the  accuracy  of  rigid 
plastic  solutions  could  be  greatly  improved  by  using  an  increased 
yield  stress  based  on  the  rate  of  generalized  strain.  However,  he 
also  concluded  that  for  the  case  of  structures  with  complex  inter- 
acting dissipation  modes,  extra  care  must  be  taken  in' assessing  the 
influence  of  the  interacting  generalized  strain  rates.  Manjoine 
(reference  8)  correlated  the  existing  experimental  data  for 
experiments  on  mild  steel  and  found  that  the  empirical  formula, 

ay  = V (1  + £/40-4>°‘20  (9) 

gave  a very  good  approximation  for  the  value  of  the  dynamic  yield 
stress. 

If  the  stationary  mode  solution,  Equation  (6) , is  substituted 
into  the  virtual  velocity  equation.  Equation  (1) , one  obtains 

W(t)  Js  Pi  (t)  4k  (x)  ds  - W(t)W(t)  Jv  m<|)i(x)cl>i(x)dv 

- Jv  Qj4jdv  , <7 8 * 10> 


Equation  (10)  can  be  rewritten  as  a non-linear,  ordinary  differential 
equation. 


where , 


P(t)W(t)  - MW ( t) W( t)  = D (W,W) 


p(t)  = Jg  Pi(t)(j)i(x)ds 


(11) 


(Ha) 


7.  Bodner,  S.  R. , "Strain  Rate  Effects  in  Dynamic  Loading  of 
Structures,"  in  Behavior  of  Materials  under  Dynamic  Loading, 
edited  by  N.  J.  Huffing ton,  Jr.,  ASME  (1965). 

8.  Manjoine,  M.  J.,  "Influence  of  Rate  of  Strain  and  Temperature  on 

Yield  Stresses  of  Mild  Steel,"  Journal  of  Applied  Mechanics, 

Vol.  II,  Trans.  ASME,  Vol.  66,  Ser.  A,  pp.  211-218  (1944). 
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M = Sv 

mcfK  (x)  $ . (x)  dv 

(Hb) 

= fv 

Q .dv 

(11c) 

The  only  step  now  remaining  to  complete  the  problem  formulation  is 
the  selection  of  an  admissible  generalized  stress  field. 

It  has  been  found  that  in  many  cases,  the  static  collapse  mode 
gives  a very  good  approximation  to  the  final  deformed  shape  of  a 
structure  loaded  dynamically  into  the  plastic  range.  Static  collapse 
modes  can  be  found  by  examining  the  displacement  pattern  for  an 
elastic-perfectly  plastic  structure  and  passing  to  the  limit  as  the 
elastic  component  of  the  displacement  becomes  negligible.  As  an 
example  of  this  process,  consider  the  simply-supported  beam  shown 
in  Figure  1.  The  elastic  deflections  are 


W ( X)  = fx  r 3 
24EI  LX 


2Lx2  + L3] 


(12) 


The  maximum  bending  moment  occurs  at  the  center  of  the  beam,  x=L/2, 
and  is  equal  to 


M(x=L/2) 


(13) 


while  the  maximum  elastic  deflection  is  equal  to 


W(x  = L/2) 


f LJ 
24EI 


(14) 


If  we  approximate  the  deflection  up  to  the  limit  point  by  the  elastic 
deflection,  the  deflection  at  the  center  under  a load  causing  a fully 
plastic  moment  at  the  center  is 


W _ 

L ~ Eh 


(15) 


For  mild  steel,  this  reduces  to 


w s 10 
L ^ h 


(16) 


Hence,  for  beams  of  appreciable  thickness,  the  elastic  deflection 
will  truly  be  negligible  in  problems  of  large  plastic  deformation. 
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As  the  load  intensity,  f,  increases  the  local  relative  deflections 
in  the  elastic  regions  will  still  remain  small,  however  large  rela- 
tive displacements  will  take  place  in  the  plastic  region  near  the 
center.  Hence,  if  we  neglect  the  curvature  changes  in  the  elastic 
regions  relative  to  the  curvature  changes  in  the  plastic  region, 
the  beam  will  deflect  in  the  form  shown  by  the  dotted  line  in 
Figure  1.  The  circle  in  the  center  represents  the  region  of  plastic 
deformation,  in  the  limit,  this  becomes  a "plastic  hinge."  Hence, 
an  appropriate  kinematically  admissible  mode  shape  for  this  problem 
would  be 


<b  fx)  = ^ 

<PW  Ll  - 2 (x-L/2)  /L, 


0£x<L/2 

L/2<x<L 


(17) 


and  a corresponding  statically  admissible  moment  field  would  be 


M 


. . TM0  (2x/L) , 0<x<L/2 

' 1 Lmq  (1  - 2 (x-L/2) /L,  l72<x<L 


(18) 


? 
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3.  ANALYTICAL  DEVELOPMENT  OF  STRUCTURAL  RESPONSE  MODEL 
3.1  Scope  of  the  Problem 

The  previous  example  of  the  limit  analysis  of  a simply- 
supported  beam  without  axial  constraints  illustrated  the  idea  of  a 
plastic  hinge  as  the  limit  of  a highly  localized  region  of  plastic 
deformation.  This  concept  will  now, be  extended  to  include  the  case 
of  bending  type  structural  elements  which  are  restrained  against 
axial  stretching.  For  such  structures,  the  dissipation  function 
will  consist  of  contributions  from  bending  at  hinge  regions  and 
from  stretching  of  cross  sections.  To  illustrate  this,  the  case 
of  a rigid-plastic  beam  with  built-in  ends  will  now  be  considered. 

For  simplicity,  the  beam  cross  section  will  initially  be 
assumed  to  be  rectangular.  Any  possible  stress  distribution  must 
be  of  the  general  form  shown  in  Figure  2,  where  z represents  the 
current  position  of  the  beam  neutral  axis.  For  a rectangular  cross 
section  of  depth  b,  the  generalized  stresses  are  an  axial  stress 
resultant 


N = / o bdz  = 2zba  (19) 

-h/2  Y y 

and  a bending  moment 

h/2  h 9 -9 

M = / CyZbdz  = Oyj  (h  - 4z  ) (20) 

In  order  to  relate  the  generalized  stresses  to  the  generalized  strain 
rates,  equations  for  the  generalized  strain  rates  as  functions  of 
the  displacement  field  and  equations  relating  the  position  of  the 
neutral  axis  in  terms  of  the  displacement  field  must  be  developed. 

To  obtain  these  relationships,  reference  must  be  made  to  a displace- 
ment mode  shape.  An  obvious  choice  for  the  displacement  mode  is 
to  use  the  static  collapse  mode.  However,  the  static  collapse  mode 
for  a beam  is  a bending  only  deformation  mode.  The  question  of 
whether  using  this  mode  is  valid  for  a problem  involving  both  bending 
and  stretching  will  be  discussed  later,  along  with  a derivation  of 
the  generalized  stress  versus  displacement  relationship  for  a beam 
with  built-in  ends, 
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N_  _ W 

N h 
Y 


(21) 

(22) 


When  Equations  (21)  and  (22) 
(or  Tresca)  yield  criterion, 


are  used  together  with  the  von-Mises 
an  interaction  relation  of  the  form, 


(23) 


is  obtained. 

Equation  (23)  is  often  also  referred  to  as  a generalized 
yield  condition.  Thus,  in  addition  to  satisfying  the  conditions 
for  being  statically  admissible,  any  generalized  stress  field  must 
also  satisfy  the  interaction  relationship.  Requiring  the  generalized 
stresses  to  follow  the  interaction  relation  is  similar  to  requiring 
the  stresses  to  follow  the  Levy-Mises  relations,  Equation  (2) , in 
classical  plasticity.  A more  detailed  discussion  of  these  points 
including  certain  restriction  on  the  generalized  strain  rates  is 
presented  by  Symonds  (reference  5) . Because  of  the  obvious 
restriction  that 


Equations  (21)  and  (22)  only  hold  for  w<h.  Beyond  this  point  one 
finds 


M 

N 


for  w > 


h 


(24) 


The  degree  to  which  interacting  generalized  stress  fields 
complicate  the  response  calculation  problem  for  even  simple  structural 
elements  is  now  obvious.  For  the  bending  only  response  calculation 
the  generalized  stress  term  is  constant  and  independent  of  displace- 
ment. The  dissipation  term  in  Equation  (11)  becomes  a linear 
homogeneous  function  of  the  generalized  strain  rate,  yielding  a linear 
differential  equation  to  be  solved.  However,  the  generalized  stress 
field  for  the  problem  involving  bending  and  stretching  is  a non-linear 
function  of  the  displacement  and  the  resulting  differential  equation 
becomes  non-linear.  The  detailed  development  of  the  kinematics  and 
dissipation  relationships  for  various  structural  members  will  now 
be  presented. 
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3.2  Beam  with  Built-in  Ends 


In  choosing  a response  mode  for  studying  the  dynamic  response 
of  structural  members  with  axial  restraints,  one  must  carefully 
assess  the  validity  of  using  a stationary  mode  to  describe  a problem 
for  which  the  generalized  stresses  and  strain  rates  are  functions 
of  the  displacement  field.  As  was  previously  mentioned,  static 
collapse  modes  have  been  found  to  give  reasonable  representations 
of  the  dynamic  collapse  mode  for  many  structures.  Sawczuk  and 
Winnicki  (reference  9)  found  that  in  reinforced  concrete  plates  under- 
going moderately  large  deflections,  the  regions  of  membrane 
deformation  were  narrow  bands  adjacent  to  the  bending  hinge  lines. 

It  was  not  until  the  deflections  became  large  that  the  overall  deflec- 
tion pattern  began  to  deviate  markedly  from  the  initial  bending 
collapse  mode.  A similar  pattern  is  found  in  examining  the  experi- 
mental results  of  Humphreys  (reference  10)  for  explosively  loaded 
beams.  The  static  collapse  mode  will  therefore  be  used  in  our 
development  of  stationary  mode  solutions  for  dynamic,  rigid-plastic, 
structural  analyses. 

The  deformation  mode  chosen  for  the  beam  is  that  shown  in 
Figure  3.  Deformations  within  the  hinge  regions  Al  and  A2  will  be 
described  by  average  quantities,  and  no  attempt  will  be  made  to 
describe  in  detail  the  stress  and  strain  variations  within  the  hinge 
regions.  The  variation  in  the  axial  stress  resultant  is  a function 
of  the  change  in  curvature  of  the  beam.  Since  we  are  assuming  a 
deformation  mode  in  which  curvature  changes  are  limited  to  very  small 
regions  of  the  beam  span,  i.e.. 


2A1  + K<  2fi> 


(25) 


it  seems  valid  to  assume  that  the  axial  stress  resultant  can  be 
accurately  approximated  as  being  constant  over  the  beam  span, 


N = N (W(t) ) (26) 


By  virtue  of  Equation  (23) , Equation  (26)  implies  that  the  average 
bending  moments  in  hinge  regions  Al  and  A2  are  equal.  Assuming 
Equation  (25)  to  remain  valid,  the  kinematics  of  the  beam  can  be 
accurately  approximated  by  neglecting  the  hinge  regions  and  assuming 
the  beam  to  behave  as  two  linear  segments.  If  S£,  represents  the 
stretching  of  One  segment, 

9~.  Sawczuk,  A.  and  Winnicki,  L.,  "Plastic  Behavior  of  Simply 
Supported  Reinforced  Concrete  Plates  at  Moderately  Large 
Deflections,"  Int.  Journal  of  Solids  and  Structures,  Vol.  1, 
pp.  97-111  (1965). 

10.  Humphreys,  J.  A.,  "Plastic  Deformation  of  Impulsively  Loaded 

Straight  Clamped  Beams,"  Journal  of  Applied  Mechanics,  Vol.  32, 
pp.  7-10  (1965) . 
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65, 


25,  + 


2 2 
W ^ VT 

5 a ^ 


25, 


(27) 


assuming  65,  <<  25,.  Differentiating  yields 


65, 


ww 

5,  * 


(28) 


Under  the  same  assumptions,  one  finds 


and 


6 


tan 


W 

5- 


^ W 
s 5, 


(29) 


(30) 


Consider  the  two  regions  Al  and  A2/2  shown  in  Figure  3.  The  two 
regions  have  identical  stress  resultants  and  angular  rotations. 

It  therefore  seems  reasonable  to  assume  they  represent  identical 
plastic  regions.  Because  of  the  assumption  of  rigid  plastic 
material,  deformation  is  limited  to  the  plastic  regions.  The  length 
of  each  of  these  plastic  regions  is  thus  65,/2.  Letting  e and  ip  be 
average  strains  and  curvature  changes  in  the  plastic  regions  yields, 


Utilizing  (27)  - (30)  yields 


£■  = 


65, 

65, 


= 


26 

6il 


e = 


* = 


25,9 

W 

45,0 

W2 


(31) 

(32) 


(33) 

(34) 


Applying  a virtual  displacement  to  the  basic  strain-curvature 
relationship  yields 


e = Z\p 


(35) 
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Combining  Equations  (33)  - (35)  one  finds 


(36) 


Since  the  position  of  the  neutral  axis  is  constrained  to  lie  in  the 
region  -h/2<z<h/2 , Equation  (36)  implies  that  the_beam  will  be  in 
a state  of  combined  bending  and  stretching  until  z = h/2,  or 
equivalently  when 


W = h 


(37) 


after  which  the  beam  responds  in  a pure  stretching  mode.  Thus,  the 
motion  of  the  beam  will  be  analyzed  in  two  phases.  The  first  phase 
covers  the  range  of  response  for  which  W _<  h and  the  beam  dissipation 
consists  of  contributions  from  both  bending  and  stretching  modes. 

When  the  displacement  equals  the  thickness  of  the  beam,  the  second 
phase,  consisting  of  pure  membrane  dissipation,  begins. 

The  relationship  between  the  position  of  the  neutral  axis 
and  the  displacement  of  the  beam  center,  Equation  (36),  is  the  key 
to  developing  the  final  form  of  the  interaction  relation.  Equations 
(21)  - (23).  Once  the  interaction  relation  has  been  found,  and  the 
various  response  modes  of  the  structure  delineated  all  that  remains 
in  the  problem  formulation  is  the  evaluation  of  the  coefficients 
(11a)  - (11c)  in  the  response  equation.  Restricting  the  loading 
to  be  a spatially  uniform  lateral  pressure,  the  generalized  force 
P(t)  is. 


P(t)  = p(t)  jj^  <j>  (x)  dsj 

where  p(t)  is  the  load  per  unit  length.  Using  Equation  (17)  as  a 
suitable  mode  form  yields 

P(t)  = p(t)  [/£  ( 2x/L)  dx  + (1  - 2 (x-&)  /L)  dxj 

or 


P(t)  = Jlp(t) 


(37) 


Following  the  same  procedure  for  M,  yields 


(38) 
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At  each  hinge  region,  the  rate  of  energy  dissipation  is 

D = /0A  [Ne  + M\|;]ds  (39) 

however,  N and  M are  treated  as  average  values  over  the  hinge  region, 
hence, 

D = N/QAeds  + M/QA^ds  (40) 


The  integral  of  the  strain  rate  over  each  hinge  region  is  just  the 
rate  of  change  of  the  length  of  the  hinge  region,  which  was  previously 
shown  to  be  equal  to  6Z/2.  The  integral  of  the  rate  of  change  of 
the  curvature  change  in  the  hinge  region  was  shown  to  be  the  rate  of 
rotation  of  the  hinge  region,  9.  Hence,  using  Equations  (28)  and 
(30)  yields  the  dissipation  rate  per  hinge  region  for  phase  1 


MWW  , MW 

N2l  + MI 


(41) 


Equation  (41)  can  be  simplified  even  further  by  using  Equations  (21) 
and  (22)  to  remove  the  implicit  dependency  of  N and  M on  W, 


or 


Di  • V 2^  + \ i - V ^ 


Di  - 


N M 

_Z Z 

2hil  h2, 


2*  M * 

W W + W . 


(42) 


Equation  (42)  represents  the  rate  of  dissipation  per  hinge  region 
for  a beam  undergoing  both  bending  and  stretching.  If  it  is 
linearized,  one  obtains  for  the  dissipation  My§  which  is  identical 
to  the  dissipation  found  for  the  bending  only  case.  To  obtain  the 
dissipation  for  the  entire  beam,  the  dissipation  per  hinge  region, 
Equation  (42)  is  multiplied  by  the  number  of  hinge  regions,  four, 
yielding 


DI  = 


2N  4M 

Z Z 

h)l  h 2l 


2* 

W W 


4M 

Z 

l 


W 


(43) 


In  phase  II,  the  entire  span  of  the  beam  is  deforming  plastically 
in  a pure  membrane  mode.  The  stretching  rate  for  each  segment  of 
the  beam  is  given  by  Equation  (28) , and  the  total  dissipation  rate 
is 
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'II 


2Ny  6 a. 


or,  equivalently. 


(44) 


'II 


2Ny»f 


(45) 


Using  Equations  (37) , (38) , (43) , and  (45)  the  governing  system  of 
differential  equations  for  the  case  of  a beam  with  built-in  ends 
is  seen  to  be: 


2 " v 

£p(t)  - |m£W  = + 

and 


Ap ( t)  - |m£W  = W,  W>h  (47) 

Equations  (46)  and  (47)  are  both  second  order  ordinary  differential 
equations  each  requiring  two  initial  conditions  to  specify  the 
solution.  For  Equation  (46) , the  initial  conditions  are 


W (o)  = 0 (48a) 

W(o)  = VQ  (48b) 

Where  VQ  is  the  assumed  magnitude  of  the  initial  modal  velocity. 

A detailed  discussion  of  how  the  modal  initial  velocity  is  chosen 
is  given  in  Appendix  A.  The  initial  conditions  for  Equation  (47) 
are  just  the  values  of  the  displacement  and  velocity  of  the  solu- 
tion to .phase  I evaluated  at  the  time  when  the  beam  begins  to  behave 
in  a pure  membrane  fashion,  i.e.,  when  W = h.  Thus,  the  complete 
solution  procedure  is  the  following: 

1.  Calculate  the  appropriate  modal  initial  conditions  for 
phase  I (Equation  46) . 

2.  Using  these  modal  initial  conditions,  generate  a solution 
Wi(t)  to  Equation  (46). 

3.  Solve  for  the  root,  t',  of  the  Equation  Wi(t’)  = 0. 

4.  Solve  for  the  root,  t*,  of  the  Equation  Wi(t*)-h  = 0. 

5.  If  t'£t*,  Wi(t')  is  the  final  deflection  of  the  beam. 

6.  If  t’>_t*,  calculate  fti(t*). 

7.  Generate  a solution  to  Equation  (47)  using  the  initial 
conditions. 


2N 

h£ 


4M_ 


h1 2 3 4 5 6 7iU 


W,  0<W<h 


(46) 
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Wl;[(o) 

= h 

(49a) 

WI;[  (o) 

= Wj(t*) . 

(49b) 

A detailed  discussion  of  the  computational  procedure  for  determining 
Wj ( t)  and  Wji  ( t)  is  given  in  Appendix  B.  An  evaluation  of  the 
accuracy  of  these  solutions  will  be  presented  in  Chapter  4. 

3.3  Beam  with  Pinned  Ends 


The  case  of  a beam  whose  ends  are  free  to  rotate  but  cannot 
translate  can  be  developed  as  a simple  extension  of  the  theory 
developed  for  the  beam  with  built-in  ends.  The  kinematic  deforma- 
tion mode  illustrated  in  Figure  3 for  the  built-in  end  case  can  be 
used  for  the  pinned  end  case  also,  however,  the  hinge  regions  denoted 
by  A]_,  at  the  two  ends  of  the  beam  now  become  zero  length  mechanical 
boundary  condition  hinges  rather  than  structural  deformation  hinge 
regions.  Following  a development  analogous  to  that  developed  for 
the  built-in  end  case,  assuming  the  normal  stress  resultant  to  be 
constant  over  the  span  implies  that  all  structural  deformation  takes 
place  in  the  central  hinge  region  until  the  entire  beam  reaches  a 
tensile  stress  state.  Because  the  assumed  displacement  mode  has 
not  changed.  Equations  (28)  and  (30)  still  hold.  However,  since 
the  stretching  SJl  is  now  concentrated  in  the  central  hinge  regions 


e = 


&JL 

61 


which  is  unchanged  from  before,  but 


* - n 


(50) 


which  implies 


• 2W 

wz 


(51) 


Since  Equation  (35)  still  holds,  one  finds 


2W  - 2W 

w - w2 


which  implies 


z = w 


(52) 
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Hence,  the  combined  bending  and  stretching  phase  of  dissipation  will 
end  when 


(53) 


and  the  motion  of  the  beam  in  the  membrane  only  phase  will  be 
identical  regardless  of  whether  the  beam  has  built-in  or  pinned  ends. 
The  dissipation  function  for  the  beam  can  now  be  stated  for  each 
phase.  Since  the  interaction  relation  is  determined  from  the 
relationship  between  the  beam  deflection  and  the  position  of  the 
neutral  axis,  the  interaction  relation  for  the  pinned  end  beam  will 
differ  from  that  found  for  the  built-in  end  beam,  Equations  (21) 
and  (22)  . Utilizing  Equation  (52)  , the  interaction  relation  for 
this  case  becomes; 


N 2W 


M_ 

My 


1-4 


(54) 

(55) 


When  Equations  (54)  and  (55)  are  combined,  a generalized  yield 
criterion  identical  to  Equation  (23)  is  obtained.  This  is  obviously 
correct  since  Equation  (2-3)  could  be  developed  solely  by  inspection 
of  Figure  2,  without  any  reference  to  the  kinematics  of  the  problem. 
The  beam  motion  can  be  described  in  two  phases  very  similar  to  that 
formulated  for  the  beam  with  built-in  ends.  However,  whereas  in 
the  previous  case  the  first  phase  of  combined  bending  and  membrane 
response  lasted  until  W = h,  for  this  case.  Equation  (54)  implies 
that  the  first  phase  will  end  and  the  beam  will  begin  to  respond 
in  a pure  membrane  manner  when  W = h/2. 


During  the  first  stage  of  the  response,  dissipation  takes 
place  in  the  central  hinge  region.  The  deformation  in  this  region 
consists  of  a stretching  of  a magnitude  261  and  a rotation  of 
magnitude  20.  Hence, 


• • • u 

Dj  = 2N5JI  + 2M0 , W<| 


(56) 


Utilizing  Equations  (54)  and  (55) , the  dissipation  function  can  be 
expressed  as 


hit 


W, 


(57) 
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Since  the  basic  kinematics  of  the  beam  are  identical  to  that  in  the 
previous  case,  the  dissipation  term  is  the  only  one  which  changes 
in  the  virtual  velocity  equation.  In  addition  upon  comparison  of 
Equations  (43)  and  (57)  , it  becomes  apparent  that  they  both  can  be 
written  in  the  general  form 

• A • • 

DI  = ax  W w + a2  W (58) 

and  thus  for  a beam  with  either  built-in  or  pinned  ends  the  response 
equation  for  phase  I can  be  written  in  the  general  form 

Jlp(t)  - |mfi,W  = a]_  W2  + a2  (59) 


Referring  again  to  the  case  of  the  beam  with  pinned  ends,  for  phase 
II  when  the  beam  is  responding  in  a pure  membrane  mode,  the  dissipa- 
tion and  hence  the  governing  differential  equation  is  identical  to 
that  of  the  beam  with  built-in  ends.  The  governing  equation  for 
phase  II  motion  of  the  beam  with  pinned  ends  is  therefore  Equation 
(47)  with  the  only  difference  being  that  motion  in  this  phase  now 
begins  when  W >_  h/2  rather  than  when  W ^ h.  If  the  phase  I to 
phase  II  transition  criterion  just  mentioned  is  used,  the  procedure 
for  calculating  the  response  of  a beam  with  pinned  ends  will  follow 
exactly  the  solution  procedure  outlined  for  the  beam  with  built-in 
ends.  The  details  of  implementing  this  solution  procedure  will  be 
developed  in  Appendix  B. 

3 . 4 Rectangular  Plates 

The  response  of  rigid  plastic  rectangular  plates  with  either 
built  in  or  pinned  ends  will  be  developed  in  a manner  similar  to 
that  used  for  beams.  The  hinge  region  dissipation  relations  found 
for  beams  will  be  shown  to  be  directly  applicable  to  plates  provided 
the  assumed  mode  shapes  satisfy  certain  criteria.  Following  the 
justification  given  for  beams,  the  static  collapse  mode  proposed 
by  Wood  (discussed  in  reference  3)  will  be  used.  This  collapse  mode 
is  composed  of  flat  sections  connected  by  hinge  regions  as  is 
illustrated  in  Figure  4.  The  solid  lines  in  the  plate  interior 
represent  hinge  regions.  The  boundary  edges  of  the  plate  represent 
either  plastic  deformation  hinge  regions  or  kinematic  boundary 
condition  hinges,  depending  on  the  boundary  conditions  applied.  For 
simplicity,  it  will  be  assumed  that  the  boundary  conditions  are  the 
same  on  all  four  sides  of  the  plate.  This  collapse  mode  has  been 
shown  by  Sawczuk  (reference  9)  to  yield  an  excellent  correlation  to 
experimental  results  on  concrete  plates.  Jones,  Uran,  and  Tekin 
(reference  11)  compared  the  same  mode  shape  to  the  permanent 


11.  Jones,  N. , Uran,  T.  0.,  and  Tekin,  S.  A.,  "The  Dynamic  Plastic 
Behavior  of  Fully  Clamped  Rectangular  Plates,"  Int.  Journal  of 
Solids  and  Structures,  Vol.  6,  pp.  1499-1512  (1970). 
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deflection  profiles  of  their  experiments  on  explosively  loaded  steel 
and  aluminum  plates  and  found  a very  good  correlation. 

If  W represents  the  displacement  of  the  center  of  the  plate, 
the  collapse  mode  for  the  quarter  of  the  plate  covering  the  region 
0<x<:B  and  0<y<B  can  be  expressed 

W = (B-y)  |, 

w ■ (L-X)  a^T' 


The  collapse  mode  for  the  other  three  quarters  of  the  plate  can  be 
obtained  from  Equation  (60)  by  symmetry  considerations.  As  can  be 
seen  from  Equation  (60) , the  collapse  mode  is  such  that  the  plate 
deflects  as  a collection  of  flat  segments  connected  along  hinge 
regions  which  themselves  can  be  represented  as  straight  lines.  But 
these  are  just  the  assumptions  invoked  in  developing  the  hinge  region 
generalized  stress  interaction  relations.  Equations  (21)  and  (22) 
or  Equations  (54)  and  (55) . Thus,  it  will  be  assumed  that  a stati- 
cally admissible  generalized  stress  field  can  be  found  which  is 
consistent  with  the  assumed  displacement  field  and  the  interaction 
relations  for  hinge  regions.  As  has  been  pointed  out  by  Jones 
(reference  11) , no  statically  admissible  stress  field  has  been  found 
for  rectangular  plates.  However,  this  does  not  mean  that  no  such 
fields  exist,  but  rather  that  no  one  has  found  a closed  form 
expression  for  such  a field.  It  will  be  assumed  that  such  a field 
does  exist  and  if  needed  could  be  found  by  numerical  methods.  This 
will  yield  a consistent  formulation  for  all  interior  hinge  regions, 
however,  it  is  not  clear  what  value  of  plate  displacement  to  associate 
with  points  on  boundary  hinges.  As  a simplifying  assumption,  the 
deflection  of  the  plate  center,  W,  will  be  associated  with  the  dissi- 
pation at  all  points  on  the  boundary.  Referring  to  Figure  4,  the 
assumed  deformation  mode,  it  is  seen  that  the  only  geometrical 
unknown  is  the  angle  <j>.  This  angle  is  assumed  equal  to  the  value 
of  it  which  gives  the  minimum  static  collapse  pressure  (reference  3) , 

(!)  = tan-1  V3  - 62  -6,e=|  (61) 


mathematically  as 


0_<x_<il^  + y tana  , 
0<y<B 


1 2_ix.iL  t 


(60) 


0<y< 


(x-^) 


—'-tan  <{> 


All  of  the  information  necessary  for  evaluating  the  terms  in 
Equation  (11)  has  now  been  developed.  From  Equation  (60) , we  find 


m 
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P(t)  = 4p(t) 


B - 


(L-x) 
tan  <l> 


i / 

L-B  tan  <|>  o 
B y tan  <j)  + {L-B  tan  4>) 

If 

o o 


B tln>]  dydx 


B 


dxdy 


(62) 


p(t)  = 


BL  B tan  $ 
2 6 


(63) 


M = 4m 


B- 


/ / 

1_L-B  tan  <\>  o 


(L-x) 
tan  4> 


L-x 


B tan  <J> 


dydx 


+ 


o o 


tan  4>  + 


(L-B  tan 


<J>) 


(64) 


M 


BL  B2  tan  4) 
3 6 


(65) 


The  dissipation  relation  will  be  developed  for  the  case  where  the 
edges  of  the  plate  are  built-in.  The  case  of  a plate  with  pinned 
edges  can  be  developed  as  a simple  extension  of  the  built-in  edge 
results  in  a manner  identical  to  that  done  for  beams.  The  dissipa- 
tion relation  will  first  be  developed  under  the  assumption  that 
W£h.  This  will  be  referred  to  as  phase  I dissipation,  as  it  was 
for  beams.  Along  a diagonal  hinge, 


and 


e = 


w 


B sin  <J> 


(66) 

(67) 
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Substituting  these  results  into  Equation  (21)  and  (22)  yields 


N.  s 


N 

2m2 

M = My  (1  - sW) 

hh22 

The  total  dissipation  along  these  diagonal  hinges  is 


\ = 4 f 


%2  " 


Ny  wV  + My  . My  wV 

l2  h h *%2 


W 


. 2 

B sxn  <j> 


ds 


D . 


Z-  N 


M. 


12  & _ ^ 


2>  W + V2 


4W 


. 2 

B sin  4> 


Following  the  same  procedure  for  the  central  hinge  yields. 


Ny  M 2 

“h  ‘ ^ ” + My 


4£,  W 
B 


and  for  the  boundary  sides  of  length  2L 


N 


M 


- W2  + 


h h2 


M 


4LW 
Y B 


plus  for  the  boundary  sides  of  length  2B 


4W 


N M 2 

DBi  = -h  - ^ W + My  tan  * 


The  total  dissipation  for  phase  I is  just 


Dt  = D , + D + Dr  + Du 

I dj  Cj  Lj  Bj 


(68) 

(69) 

(70) 

(71) 

(72) 

(73) 

(74) 

(75) 


However,  by  inspection,  it  is  obvious  that  Equation  (75)  can  be 
written  in  a form  identical  to  Equation  (58) . Hence,  the  solution 
procedure  for  phase  I motion  of  a rectangular  plate  is  identical 
to  the  solution  procedure  outlined  in  Appendix  B.  When  the  center 
deflection  of  the  plate  exceeds  the  plate  thickness,  phase  II  of 
the  plate  motion  begins.  If  s*  is  defined  as  the  distance  along 
the  diagonal  hinge  at  which  the  deflection  is  equal  to  h,  then 
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and 


D 


II 


s*  = 


hi, 

£ 

w 


= 4 


'' 

/ " 


y h 


W + - M, 


2 2 
S W 


w 


y y h2Jl2  B sin2 


ds 


(76) 


f n — i 

J y l2  B sin2  4> 


ds 


(77) 


which  gives 


‘II 


— 

r 0 

- 

N 

2M  h&0  N hi, 

y 2 w + 

y 2 y 2 

1 

2 W + 

3 6 

W 

- 

_ 

4W 


B sin  <}> 


(78) 


Up  to  this  point,  all  of  the  dissipation  relations  have  been  written 
in  terms  of  the  generalized  stresses,  Ny  and  My.  However,  the  inter- 
action relations  were  developed  with  the  inherent  assumption  that 
the  cross  section  was  homogeneous  and  uniform.  Under  these  conditions 

°vh2 

Ny  = CTyh  and  My  = 


If  these  relations  are  substituted  into 


Equation  (78) , the  term  multiplying  W-1  becomes  zero  and  the 
dissipation  for  the  diagonal  hinges  becomes 


2N  Z7m 
x - 

B sin2<J) 


(79) 


The  total  plate  dissipation  is 


2Sj2 

B sin2  c|> 


+ 


4 

tan  <S> 


NyWW 


(80) 


and  the  resulting  differential  equation  is  of  a form  identical  to 
Equation  (47) , so  that  the  solution  procedure  for  plates  is  identical 
to  that  outlined  for  beams. 


In  summary  it  can  be  seen  that  the  formulations  for  rectan- 
gular cross  section  beams  and  solid,  homogeneous  rectangular  plates 
all  yield  identical  systems  of  differential  equations  to  be  solved. 
This  system  can  be  written  in  the  general  form 
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ll 

Phase  I,  W<h  (built-in  ends)  or  (pinned  ends) 

A1  p(t)  - A2  W = A3W2  + A4  (81) 

Phase  II,  W>h  (built-in  ends)  or  W>^  (pinned  ends) 

A1  p(t)  - A2  W = A5W  (82) 

In  appendix  B a general  solution  procedure  following  the  form  out- 
lined for  a beam  with  built-in  ends  is  developed  in  detail. 

3 . 5 Extension  of  Method  to  General  Compact  Sections 

The  previous  development  for  beams  was  limited  to  beams  of 
rectangular  cross  section  because  of  the_form  of  the  interaction 
relation  used.  If,  general  forms  for  N(z)  and  M(z)  can  be  found 
for  a given  compact  section,  the  relations  already  developed  for 
z (W)  and  the  generalized  strain  rates  could  then  be  used  to  complete 
the  formulation.  However,  in  general  the  resulting  system  of 
differential  equations  will  not  be  of  the  general  form  developed 
in  the  previous  cases.  How  such  a procedure  is  implemented  will  be 
demonstrated  by  developing  the  governing  equations  for  the  I-beam 
section  shown  in  Figure  5.  Let  R represent  the  geometrical  region 
covered  by  the  beam  cross  section  in  the  x-z  plane.  Then 

t 

N = JJ  crdxdz  (83) 

R 

M = //  azdxdz  (84) 

R 

Let  Ri  be  the  region  R (z<z)  and  R2  be  the  region  R (z>z) . By  use 
of  Green's  theorem  (reference  12),  Equations  (83)  and  (84)  can  be 
written 

N = ^3R  axdz 
f 

M = J 3-  axzdz 


12.  Courant,  R.  and  Hilbert,  D. , Methods  of  Mathematical  Physics, 
Volume  1,  Interscience  Publishers,  Inc.,  (1966). 
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or  equivalently 


N ■ l 

* ■ °y L 


xdz  - a. 


xzdz  - a. 


xzdz 


Assuming  R can  be  decomposed  into  a series  of  n linear  segments  of 
the  form 


xi  = ai  + b.  z 


N and  M can  then  be  expressed  as  functions  of  z 


N ill?  [ai(Z2i  zli)  + 2 (z2i2  zli'L)] 

i=j+l  °Y  [ai(Z2i  " Zli)  + (Z2i  " Zli  {8 

3 fai  2 2 bf  3 3-1 

M = ^ ay  [-f  (z2.2  - z,.2)  + (z2i  - zn3)] 

V rai  / 2 2,  , bi  , 3 3,1  .. 

i=j+3.  °y  L 2 Z2i  zli  ] 3 (z2i  zli  (* 

where  we  use  a numbering  convention  for  the  boundary  segments  such 
that 


Z11  = z2j  “ zl(j+l)  = z2n  ~ z ’ 

The  preceeding  development  holds  for  any  region  R.  When 
Equations  (88)  and  (89)  are  evaluated  for  the  I-beam  of  Figure  5, 
the  following  relations  are  found. 

Region  1,  0<z<H 


N = ClZNy 
where,  CL  = ^ 
M = My  - C2z2My 


bHT-H(b-d) 
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where,  C 2 


d 

b(HT2-H2) +dH2 


(92a) 


Region  2,  H^ZjcHT 


= H(b-d)-bz 
H (b-d) -bHT  y 


(93) 


M = 


2 —2 
bHT  - bz 


b(HT2-H2)  + dH 


2 My 


(94) 


Utilizing  the  generalized  strain-displacement  relations  developed 
for  beams  with  built-in  ends,  Equations  (91) -(94)  become  interaction 
relations  for  the  generalized  stresses  versus  beam  displacement. 

When  W>2HT,  the  beam  will  respond  in  a pure  membrane  mode  with 
generalized  stresses,  N = Ny  and  M = 0.  The  resulting  system  of 
differential  equations  is  obtained. 


Phase  I,  Region  1,  0<w<2h 


2 ••  W2  ^ 

Jlp(t)  - fmilW  = (C^  - C2My)  \ + -j 


Phase  I,  Region  2,  2H<W<2HT 


Nyb 


Jlp(t)  - 3m£W  - bHT_H(b-d) 


Mvb  w2 

b(HT2-H2)+dH2  1 


(95) 


2H(b-d)Ny  w 
bHT-H(b-d)  i 

4bHT2M 

(b (HT2-H2) +dH2) H 


(96) 


Phase  II,  W>2HT 


ilp  ( t)  - |milW  = N W (97) 

^ y 

This  system  is  slightly  different  than  the  general  system  found  to 
be  applicable  for  both  rectangular  cross-section  beams  and  solid, 
homogeneous,  rectangular  plates.  Equations  (95)  and  (97)  form  a 
system  identical  in  form  to  the  general  system  previously  found. 

In  addition  to  adding  a stage  of  beam  response  which  must  be  included 
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in  the  solution  procedure.  Equation  (96)  is  of  a more  complex  nature 
than  the  other  equations  which  have  been  developed.  As  will  be  seen, 
the  solution  algorithm  developed  in  Appendix  B is  general  enough 
to  cover  the  system  of  equations  for  the  I-beam,  Equations  (95) , 

(96) , and  (97) . 
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4.  SAMPLE  RESULTS  AND  COMMENTS  ON  STRUCTURAL  RESPONSE  MODEL 

In  this  chapter,  a series  of  sample  problems  solved  using  the 
methods  developed  in  the  previous  chapters  will  be  presented.  The 
purpose  of  this  chapter  is  twofold,  first  to  illustrate  the  accuracy 
of  the  solution  techniques  which  have  been  developed  and  second 
to  point  out  several  significant  points  which  designers  of  suppressive 
shields  must  take  into  account.  For  simplicity,  in  discussing  the 
results  of  calculations,  the  methods  developed  in  this  report  will 
be  referred  to  as  the  virtual  velocity  method.  This  is  not  meant 
to  be  a completely  descriptive  title,  but  rather  just  a convenience 
for  specifying  the  methods. 

4 . 1 Impulsively  Loaded  Rectangular  Cross-Section  Beam 

It  is  now  appropriate  to  establish  the  accuracy  of  the 
virtual  velocity  method  when  compared  to  published  experimental 
results.  There  is  quite  an  amount  of  data  available  for  the  rigid 
plastic  response  of  beams  without  axial  constraint  to  impulsive 
loads.  However,  very  little  data  is  available  for  the  case  of  beams 
with  axial  constraint. 

Consider  as  an  example,  a rigid  plastic  beam  with  built- 
in  ends  (with  axial  constraint)  subjected  to  an  impulsive  load. 

The  virtual  velocity  method  as  described  in  Chapter  3 was  used  to 
predict  the  permanent  deformation  in  the  beam.  The  results  of  this 
analysis  are  compared  with  the  experimental  data  of  Humphreys 
(reference  10)  in  Figure  6. 

The  virtual  velocity  method  agrees  with  the  experimental 
data  for  mild  steel  quite  well  throughout  the  range  of  deflections. 
However,  the  agreement  with  the  data  for  the  4130  steel  is  not  as 
satisfactory.  This  occurs  because  the  4130  steel,  a chromium  moly- 
bdenum steel,  does  not  exhibit  the  same  variation  of  material 
properties  withstrain  rate  as  mild  steel. 

The  virtual  velocity  method  appears  to  agree  with  the 
experimental  data  better  than  the  computationally  more  complex  method 
of  Symonds  and  Mentel  (reference  13)  for  the  region  of  membrane 
only  response.  One  explanation  for  this  is  the  inclusion  of  strain 


13.  Symonds,  P.  S.  and  Mentel,  T.  S.,  "Impulsive  Loading  Of  Plastic 
Beams  with  Axial  Constraints,"  Journal  Of  The  Mechanics  and 
Physics  of  Solids,  Vol.  6,  pp.  106-202  (1958). 
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rate  effects  which  Symonds  and  Mentel  ignored.  For  the  region  of 
small  deflections  i.e.  S/h  £ 1/2  the  bending  only  solution  yields 
acceptable  results.  This  is  reasonable  because  in  that  region  the 
response  is  mainly  due  to  bending.  For  larger  displacement,  i.e. 

S/h  £ 2,  as  the  membrane  response  becomes  dominent,  the  bending  only 
solution  gives  totally  erroneous  results. 

4 . 2 Comparison  between  I-Beam  Section  and  Rectangular  Section 
Beam  of  Equal  Mass  Per  Unit  Length 

Designers,  choosing  beam  sections,  use  structural  shapes 
which  are  more  complex  than  rectangular  shapes.  For  example,  in 
the  quarter-scale  category  1 suppressive  shield,  standard  section 
I-beams  are  used  to  construct  the  wall.  For  most  applications,  a 
compact  section  such  as  a standard  section  I-beam  gives  a more 
efficient  design  than  does  a rectangular  section.  Since  many  simpli- 
fied design  methods  do  not  accurately  model  the  membrane  - bend- 
ing interaction  in  their  treatment  of  I-beams,  it  is  important  to 
determine  the  difference  in  the  response  characteristics  of  an  I-beam 
section  and  a rectangular  cross-section  beam  of  equal  mass  per  unit 
length.  Since  both  beams  are  of  the  same  density,  this  means  that 
both  beams  have  equal  cross  sectional  area.  In  addition,  the  width 
of  the  rectangular  section  was  chosen  to  be  equal  to  the  width  of 
the  I-beam  so  that  equal  pressure  would  cause  equal  loads  on  the 
two  beams.  For  this  comparison  the  depth  of  the  I-beam  is  consi- 
derably greater  than  that  of  the  rectangular  section  beam.  Figure  7 
shows  the  results  obtained  from  such  a virtual  velocity  method 
analysis.  Notice  that  there  is  considerably  less  deflection  in  the 
I-beam  than  in  the  rectangular  cross-section  beam  for  any  given 
impulse. 

Since  the  I-beam  is  much  deeper  than  the  rectangular  beam 
its  bending  rigidity  is  significantly  greater  than  that  of  the  rectan- 
gular beam.  Consequently,  bending  controls  the  response  of  the  I- 
beam  for  a wider  range  of  deflections.  As  a result  the  I-beam 
deflects  less  for  any  given  impulse  than  does  the  rectangular  beam. 

Once  both  beams  reach  the  pure  membrane  realm,  their  response 
characteristics  are  similar.  This  may  not  be  readily  apparent;  how- 
ever, recall  that  the  governing  equation  of  motion  for  the  pure 
membrane  response  is  controlled  by  the  cross  sectional  area  and  the 
yield  stress.  Also  recollect  that  for  this  illustrative  case  both 
beams  have  equal  cross  sectional  area  and  yield  stress.  Therefore, 
it  is  evident  that  in  the  pure  membrane  realm  the  governing  equation 
of  motion  for  both  beams  are  identical  and  the  solution  to  such 
equations  should  yield  similar  results. 

4 . 3 I-Beam  and  Rectangular  Beam  of  Equal  Properties 

It  is  important  to  determine  the  effect  of  the  bending  - 
membrane  interaction  on  the  plastic  response  of  various  structural 
shapes  with  identical  plastic  properties.  That  is,  because  different 
structural  shapes  have  different  interaction  relationships  the 
structural  response  of  each  shall  differ  accordingly.  Therefore  a 
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comparative  study  was  performed  between  a rectangular  cross-section 
beam  and  I-beam  section  with  identical  plastic  yield  moment  (M  ) 
and  stress  resultant  (N^) . ^ 

In  order  to  obtain  section  properties  identical  to  a standard 
shape  I-beam  (S  3 x 5.7) , the  rectangular  beam  must  be  some  2.5  times 
as  deep  as  the  I-beam  and  only  15%  as  wide.  Figure  8 shows  the 
results  of  a comparison  between  a rectangular  section  and  a standard 
shape  I-beam  loaded  impulsively  into  the  plastic  range. 

The  responses  of  both  are  nearly  identical  up  to  a deflec- 
tion of  0.05  m (2  in).  However,  beyond  that  point  the  response  curves 
begin  to  differ.  The  difference  icnreases  until  a deflection  equal 
to  the  thickness  of  the  B-beam  is  reached.  At  that  point  the  deflec- 
tion in  the  I-beam  is  some  10%  greater  than  that  of  the  rectangular 
section. 


From  this  fact,  one  would  imagine  that  the  rectangular 
section  beam  would  be  the  more  efficient  section.  However,  because 
the  rectangular  section  is  very  high  and  narrow,  over  six  rectangular 
section  beams  would  be  required  to  occupy  the  same  space  as  one  I- 
beam.  Consequently  the  weight  of  the  structure  would  increase  by 
that  amount,  not  an  efficient  use  of  steel. 

As  the  deflection  increases  to  beyond  0.1m  (5  in)  the 
difference  between  the  two  sections  remains  nearly  constant  at  about 
0.36m  (14  in).  Since  both  sections  have  identical  cross  sectional 
area  and  yield  stress  the  similarity  in  response  is  expected.  This 
agrees  with  the  previous  example. 

As  can  be  seen  for  larger  deflections,  not  including  the 
bending-membrane  interaction  underpredicts  the  final  deflection. 

In  the  case  of  large  plastic  deformation  this  difference  may  be 
quite  significant.  For  design  purposes  that  effect  of  underpredicting 
the  deflection  is  not  conservative.  Therefore,  the  bending-membrane 
interaction  is  significant  to  the  overall  structural  response  and 
should  be  included  in  the  response  formulation. 

In  a similar  manner  a comparison  can  be  made  between  a wide 
flange  I-beam  and  a rectangular  section  beam  with  identical  section 
properties.  In  order  to  obtain  section  properties  identical  to  the 
wide  flange  I-beam  (W  10  x 77) , the  rectnagular  beam  must  be  some 
70%  deeper  and  only  15%  as  wide. 
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The  trend  of  the  results  is  similar  to  those  found  for  the 
standard  shape  I-beam.  The  response  of  both,  shown  in  Figure  9, 
is  nearly  identical  up  to  a point.  However,  beyond  that  point  the 
response  begins  to  differ.  The  difference  increases  until  a deflec- 
tion of  about  0.41m  (16  in)  is  reached.  At  that  point  the  deflection 
of  the  wide  flange  is  about  10%  greater  than  that  of  the  rectangular 
beam.  The  difference  then  remains  fairly  constant  beyond  that  point 
for  the  reason  discussed  previously. 

Again  it  must  be  emphasized  that  for  larger  deflections, 
not  including  the  effects  of  the  bending-membrane  interaction  under- 
predicts the  deflection,  a result  that  may  be  quite  significant  since 
it  is  not  conservative  design  practice.  Likewise  this  is  not  to 
imply  that  the  rectangular  section  is  a more  efficient  structural 
element  because,  as  before,  over  six  would  be  required  to  occupy 
the  same  space.  Therefore  the  structural  weight  would  increase  by 
that  amount. 

4 . 4 Flat  Plates,  Impulsively  loaded 

The  virtual  velocity  method  was  applied  to  the  problem  of 
a rigid  plastic  rectangular  flat  plate  with  built-in  ends  which  was 
loaded  impulsively.  The  results  of  the  virtual  velocity  method  are 
compared  with  the  experimental  results  of  Jones,  Uran  and  Tekin 
(reference  11)  in  Figure  10.  There  is  excellent  agreement  with  the 
experimental  results  for  all  ranges  of  deflections. 

4.5  Finite  Duration  Load 


There  is  an  absence  of  experimental  data  on  the  plastic 
response  of  beams  and  plates  subjected  to  finite  duration  loads. 
However,  in  the  case  of  an  accidental  explosion  within  a suppressive 
shield,  the  loads  will  be  caused  by  a shock  load  as  well  as  a quasi- 
static pressure  load.  It  is  important  to  determine  the  effect  of 
venting  on  the  structural  response  of  suppressive  shield  structures. 
To  this  end,  analytical  studies  utilizing  the  virtual  velocity  method 
were  performed  to  determine  the  net  effect  of  venting. 

Figures  11  and  12  shows  the  loading  pulses  for  the  vented 
and  unvented  quarter-scale  category  1 suppressive  shields 
respectively.  For  the  vented  condition  the  gas  pressure  is  assumed 
to  vent  from  1.38  MPa  (200  psi)  to  atmospheric  in  about  50  ms,  while 
for  the  unvented  condition  the  gas  pressure  remains  constant  at 
1.38  MPa  (200  psi) . 

A comparison  of  the  permanent  deformation  for  the  vented 
and  unvented  cases  is  made  in  the  tables  included  with  Figures  11 
and  12.  As  can  be  seen,  the  net  effect  of  venting  on  the  overall 
structural  response  is  slight.  The  effect  of  venting  would  be  more 
pronounced  if  the  vent  time  was  decreased.  However  if  this  was 
done  the  transmitted  shock  would  be  much  more  severe. 
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4.6  Quarter-Scale  Category  1 Suppressive  Shield  with  One  Ring 


The  loading  pulse  used  for  the  preceeding  venting  calculation 
was  only  approximate.  As  described  in  Chapter  5,  the  hydrocode  TUTTI 
was  used  to  predict  the  loading  function  on  the  wall.  The  loading 
function,  output  by  TUTTI,  was  used  as  input  for  the  virtual  velocity 
method  to  predict  the  plastic  response  of  the  structure. 

A hydrocode  calculation  was  performed  for  the  detonation 
of  a 22  kg  (48  lb)  charge  of  high  explosive  at  the  geometric  center 
of  a quarter-scale  category  1 suppressive  shield.  The  resulting 
loading  function  is  shown  in  Figure  13.  The  inside  beam  sees  the 
initial  shock  load  which  then  attenuates.  The  outside  beam  does 
not  see  the  initial  shock,  however,  as  the  gas  vents  through  the 
beams  the  load  rises  and  then  decays.  At  the  end  of  the  pulse, 
the  decay  of  the  outside  pulse  is  much  less  than  that  of  the  inside. 

The  response  of  the  I-beams  to  the  loading  pulse  described 
above  is  shown  in  Figure  14.  As  expected,  the  inside  beam  deflects 
more  than  the  outside  one.  This' occurs  because  the  initial  shock 
drives  the  inside  beam  while  the  venting  gas  pressure  drives  the 
outside  beam.  The  velocity  record,  given  in  Figure  15,  shows  the 
same  trend.  That  is,  the  velocity  of  the  inside  beam  rises  very 
rapidly  then  decays  while  the  velocity  of  the  outside  beam  rises 
less  rapidly  and  decays  less  rapidly. 

The  fact  that  the  inside  beam  comes  to  rest  before  the 
outside  one  can  be  explained  by  examining  the  loading  function, 

Figure  13.  As  can  be  seen,  for  later  time  in  the  load  curves, 

(t  _>  1.0  ms)  , the  load  on  the  outside  beam  decays  much  less  rapidly 
than  that  on  the  inside  beam.  This  occurs  because  the  flow  of  the 
gas  is  choked  at  the  outside  beam.  The  pressure  drop  therefore  is 
constant  and  the  load  is  controlled  by  the  upstream  flow.  The  effect 
is  to  cause  the  load  on  the  outside  beam  to  be  nearly  constant  in 
time.  This  effect  is  seen  in  the  response,  by  the  fact  that  the 
inside  beam  comes  to  rest  before  the  outside  one. 

4 . 7 The  Effect  of  Multiple  Shocks 

It  is  impossible  to  evaluate  the  effect  of  timing  of  the 
load  application  in  some  approximate  methods  of  predicting  plastic 
response.  However,  the  virtual  velocity  method  not  only  enables 
the  transient  response  to  be  determined,  but  also  enables  the 
influence  of  multiple  shocks  to  be  determined.  , 

As  an  example,  consider  the  loading  function  shown  in 
Figure  16.  This  represents  a possible  loading  function  in  which 
multiple  shocks  are  included.  Notice  that  a second  shock  arrives 
about  2.3  ms  after  the  first  shock.  The  velocity  response  which 
is  shown  in  Figure  17  is  particularly  interesting.  The  first  part 
of  the  record  is  similar  to  the  previous  velocity  record,  Figure  15. 
That  is  the  velocity  rises  rapidly  then  decays.  However,  when  the 
second  pulse  arrives,  the  velocity  again  rises  rapidly  then  decays. 
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The  inside  beam  had  nearly  stopped  when  the  second  loading  spike 
arrived.  It  then  rose  to  a higher  velocity  then  the  outside  beam. 

As  before  the  rata  of  decay  of  the  load  for  the  outside  beam  is  less 
than  that  for  the  inside,  and  the  outside  beam  again  comes  to  rest 
after  the  inside  one. 

The  deflection,  shown  in  Figure  18,  exhibits  the  same 
characteristics  as  the  velocity  record.  That  is,  when  the  second 
pulse  arrives  the  beams  deflect  more  rapidly.  From  the  velocity 
record,  this  is  what  would  be  expected. 

For  this  example,  the  second  pulse  is  only  about  one  third 
the  magnitude  of  the  first.  If  the  second  pulse  was  of  equal 
strength  to  the  first  the  effect  on  the  response  would  be  much  more 
significant.  Thus  it  can  be  seen  that  the  effect  of  multiple  wall 
reflections  is  one  which  must  be  included  in  the  analysis  if  it 
exists. 

4 . 8 Comments 

It  has  been  demonstrated  that  stationary  mode  solutions 
developed  from  the  principle  of  virtual  velocities  and  including 
the  effect  of  finite  deflections  and  dependence  of  yield  stress  on 
strain  rate  offer  the  designer  of  suppressive  shield  structures  a 
powerful  design  tool.  Through  the  use  of  solutions  generated  by 
such  procedures,  it  has  been  shown  that  conservative  design  proce- 
dures should  in  general  include  the  actual  form  of  the  interaction 
relation  between  the  generalized  stresses  and  the  time  variation 
of  the  loading  function. 
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5.  HYDROCODE  COMPUTATION  OF  SHOCK  ATTENUATION  FOR  SEVERAL 

CONFIGURATIONS  OF  I-BEAMS  AND  PERFORATED  PLATES 

5.1  Shock  Attenuation  By  A Single  Perforated  Plate 

In  this  chapter,  a method  is  developed  for  calculating  the 
details  of  the  flow  of  a shock  wave  through  various  configurations 
of  I-beams  and  perforated  plates. 

Before  considering  real  shield  configurations,  we  will 
examine  the  simple  case  of  a step  shock  hitting  a perforated  plate 
head-on.  This  case  will  be  used  as  a check  on  the  calculations. 

5.1.1  Simple  Theory 

A simple  theory  for  a single  plate  is  easily  derived.  For 
a small  vent  fraction,  we  can  assume  that  the  shock  relection  is 
identical  to  that  for  a rigid  plate.  The  reflected  shock  forms  the 
reservoir  for  flow  thru  the  orifices.  The  jets  of  gas  are  assumed  to 
mix  with  the  downstream  gas  and  conservation  of  momentum  gives  a 
downstream  particle  velocity  from  which  the  other  downstream  shock 
properties  may  be  obtained. 

5.1.2  TUTTI  Computer  Code 

TUTTI  is  a two-dimensional  Eulerian  finite-difference 
computer  program  written  at  the  Naval  Surface  Weapons  Center,  White 
Oak  Laboratory  (NSWC/WOL)  based  on  the  FLIC  (Fluid-in-Cell)  method 
of  Gentry,  et  al  (reference  14) . FLIC  originally  had  one  material; 
TUTTI  has  two.  We  have  extended  TUTTI  to  have  venting  walls  along 
any  cell  boundaries.  The  flow  thru  the  wall  is  calculated  from  the 
conditions  across  the  wall,  with  the  higher-pressure  side  being  the 
reservoir  for  flow  thru  the  orifices.  The  flow  equations  are  for 
subsonic  or  sonic  flow,  as  appropriate  (reference  15) . These  venting 
walls  make  practical  the  calculation  of  flows  thru  perforated  plates; 
most  perforated  plates  have  holes  that  are  small  compared  to  the 
distance  between  holes  and,  for  multiple  plates,  the  distance 
between  plates.  The  computer  storage  and  running  times  would  be  too 
large  if  the  individual  holes  were  modelled.  Structures  like  I-beams 


14.  Gentry,  R.  A.,  Martin,  R.  E.  and  Daly,  B.  J.,  "An  Eulerian 
Differencing  Method  for  Unsteady  Compressible  Flow  Problems," 
Journal  of  Computational  Physics,  Vol.  1,  pp. 87-118  (1966). 

15.  Proctor,  J.  F.,  "Internal  Blast  Damage  Mechanics  Computer 
Program,"  Naval  Ordnance  Laboratory  NOLTR  72-231,  AD759002, 
Silver  Spring,  Maryland  (Aug  1972) . 
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can  be  modelled  directly  with  solid  wall  segments  because  they  do  not 
mix  small-  and  large-scale  parts.  Uniformly  vented  infinite  plates 
normal  to  the  flow  reduce  to  a one-dimensional  problem  solved  with  a 
single  r(ow  of  cells  in  TUTTI . 

There  are  two  ways  to  put  in  the  reflection  of  the 
incident  shock:  1)  Put  in  the  incident  shock  before  reflection  and 

calculate  the  reflection  as  part  of  the  problem,  and  2)  Make  a 
separate  calculation  of  the  shock  reflection  from  a rigid  wall  and 
use  these  pressures,  densities,  and  internal  energies  as  time- 
dependent  inflow  in  the  TUTTI  calculation.  Method  1)  is  more  rigorous 
and  must  be  used  when  the  vent  fraction  is  large.  Method  2)  is  usable 
when  the  vent  fraction  is  small  (10  percent  or  less)  and  the  reflected 
shock  from  the  perforated  plate  is  nearly  the  same  as  that  from  a 
rigid  plate. 

5.1.3  Experimental  Data 

Dain  and  Hodgson  (reference  16)  did  experiments  with  a 
C02~filled  shock  tube  in  which  they  used  a perforated  plate  to  produce 
weak  shocks.  The  simple  theory  described  above  fits  their  data  quite 
well  (Figure  19) , as  does  a calculation  with  TUTTI,  shown  as  a 
triangle.  They  also  derived  a simple  theory  (not  shown  here)  with 
different  assumptions;  it  gives  somewhat  higher  transmitted  pressure 
than  does  our  theory. 

Kingery  and  Coulter  (reference  17)  did  shock  tube 
experiments  with  single  perforated  plates.  Figure  20  shows  their 
results  for  10  percent  venting.  TUTTI  results  are  in  reasonable 
agreement  with  the  data.  The  simple  theory  gives  transmitted 
pressures  that  are  too  low  for  strong  incident  shocks.  Figure  21 
shows  the  experimental  and  calculated  transmissions  as  a function  of 
vent  fraction.  The  agreement  between  experiment  and  TUTTI 
calculations  is  good  for  vent  fractions  below  10  percent. 

The  conclusion  from  these  comparisons  with  experiment  is 
that  for  the  low  vent  fraction  range  of  interest  in  suppressive 
shielding,  TUTTI  gives  acceptable  accuracy.  Although  the  simple 
theory  might  be  modified  to  give  better  accuracy  than  it  now  has  for 
strong  shocks,  it  cannot  give  the  transient  loads  on  the  structure 
which  are  part  of  the  suppressive  shielding  problem. 

5 . 2 The  BRL  4-Plate  Shock  Tube  Experiment 

Calculations  were  made  of  the  recent  BRL  0.102  m (4.00  in) 
diameter  shock  tube  experiment  (reference  18)  with  four  plates,  each 


16.  Dain,  C.  C.  and  Hodgson,  J.  P.,  "Generation  of  Weak  Shock  Waves 
in  a Shock  Tube,"  Aero.  Quart.,  pp.  101-108,  (May  1974). 

17.  Kingery,  C. , and  Coulter,  G. , "Airblast  Attenuation  by  Perforated 
Plates,"  Ballistic  Research  Laboratories  Interim  Memorandum  Report 
No.  338  (Feb  1975). 

18.  Kingery,  C. , Ballistics  Research  Laboratories,  to  be  published. 
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with  five  0.0127  m (0.500  in)  holes  and  a vent  factor  of  0.078  per 
plate.  The  plate  thickness  and  spacing  were  both  0.00635  m 
(0.25  in).  The  incident  shock  was  0.483  MPa  (70  psi)  and  the 
reflected  peak  pressure  on  the  first  plate  was  2.14  MPa  (310  psi). 

We  took  the  measured  pressure  vs.  time  data  on  the  first  plate  as 
input  and  calculated  the  pressures  on  the  other  plates  with  TUTTI; 
the  results  are  shown  in  Figure  22.  The  calculated  and  experimental 
results  are  compared  in  Figure  23.  The  calculated  pressures  on  plates 
2,  3,  and  4 are  larger  than  the  experimental  values  but  the  transmitted 
pressure  is  smaller  than  measured. 

5 . 3 Shock  Attenuation  By  The  0-4  1/16-Scale  Shield 

The  0-4  shield  consists  of  four  identical  steel  plates 
perforated  with  a large  number  of  small  round  holes  occupying  ten 
percent  of  each  plate  area.  The  plate  spacing  is  0.0127  m (0.500  in). 

5.3.1  Limitations  of  Simple  Theory 

The  simple  theory  described  above  could  be  tried  for 
multiple  plates  by  letting  the  downstream  shock  become  the  incident 
shock  for  the  next  plate.  However,  this  would  give  only  the  very 
first  shock  to  come  out  of  the  last  plate  and  would  not  account  for 
the  filling  of  the  interplate  volumes  or  the  multiple  reflections 
between  plates  and  their  downstream  coalescence.  The  final  steady- 
state  flow  produced  by  an  incident  step  shock  is  driven  by  the 
steady-state  pressure  in  the  last  volume  and  does  not  correspond  to 
the  "first-thru"  shock. 

A steady-state  flow  theory  for  multiple  plates  could  be 
developed.  This  would  involve  solving  for  the  interplate  and  down- 
stream pressures  with  the  conditions  that  the  mass  flux  be  equal  thru 
all  plates  and  consistent  across  each  plate  with  the  subsonic/sonic 
flow  equations.  The  density  vs.  pressure  relation  would  be  obtained 
from  isentropic  expansion  from  the  first  reservoir  (the  first 
reflected  shock  conditions) . The  effect  of  heating  by  the  multiple 
shocks  between  plates  disappears  because  this  gas  is  swept  downstream 
and  replaced  by  isentropically  expanding  flow.  This  direction  was 
not  pursued  here  because  the  actual  problems  we  are  interested  in  are 
not  for  step  shocks,  but  for  decaying  shocks  from  explosions. 

5.3.2  Calculations  for  3.1  MPa  Step  Shock 

Figures  24  and  25  show  the  interaction  of  a 3.1  MPa 
incident  step  shock  with  the  0-4  shield.  Figure  24  shows  the 
pressures  on  the  plates  vs.  time;  the  oscillations  of  the  uppermost 
curve  are  artifacts  of  the  calculation.  The  oscillations  of  the  other 
curves  are  partly  artifacts  and  partly  the  ringing  of  the  first  shock 
between  adjacent  plates.  Note  the  slow  rise  times:  the  peak  of  a 
quickly-decaying  incident  shock  would  be  strongly  clipped  by  this 
shield.  The  larger  the  space  between  plates,  the  longer  the  fill 
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time  and  the  greater  the  attenuation  of  the  peak.  Figure  25  shows 
the  pressure  distribution  at  0.2  ms  after  impact.  The  transmitted 
wave  has  not  shocked  up  yet. 

5.3.3  Calculations  for  3.1  MPa  Decaying  Shock 

Figures  26,  27  and  28  show  the  interaction  of  the  decaying 
shock  from  a 0.22  kg  pentolite  charge  detonated  in  the  center  of  a 
cubical  0-4  shield.  The  plate  1 upstream  loading  conditions  were 
obtained  from  a separate  one-dimensional  calculation  of  a charge  in 
a rigid  sphere.  The  back-reflection  of  the  reflected  shock  from  the 
HE  product  gas-air  contact  surface  is  included.  Figure  26  shows  the 
pressure  vs.  time  on  each  plate.  Note  that  the  fill  times  of  the 
chambers  between  plates  are  large  compared  to  the  decay  time  of  the 
loading  shock;  this  gives  considerable  attenuation  of  the  peak 
pressure.  Figure  27  shows  that  the  pressures  in  the  chambers  are 
practically  equal  and  that  the  last  plate  bears  the  greatest  load. 
Figure  28  shows  the  calculated  peak  pressure  just  outside  the  shield; 
it  more  or  less  agrees  with  the  experimental  data. 

5 . 4 Shock  Attenuation  By  The  T-5  1/16-Scale  Shield 

The  T-5  shield  consists  of  interleaved  I-beams.  In  the 
1/16-scale  tests,  the  widths  of  the  four  slits  through  which  the  flow 
passes  are  0.011,  0.0033,  0.0033,  and  0.011  m;  the  corresponding 
vent  fractions  are  0.156,  0.094,  0.094,  and  0.156.  Figure  29  shows 
the  zoning  for  the  calculation  of  flow  through  this  shield. 

5.4.1  Calculations  for  3.1  MPa  Step  Shock 

Figures  30  and  31  show  the  interaction  of  a 3.1  MPa 
incident  step  shock  with  the  T-5  shield.  Figure  30  shows  the 
pressures  in  the  chambers  within  the  shield  vs.  time.  Figure  31  shows 
the  pressure  distribution  within  the  shield  at  times  of  0.153  and 
0.311  ms  after  shock  impact. 

5.4.2  Calculations  for  3.1  MPa  Decaying  Shock 

Figures  32,  33,  and  34  show  the  interaction  of  the  decaying 
shock  from  a 0.22  kg  pentolite  charge  detonated  in  the  center  of  a 
cubical  T-5  shield.  The  upstream  loading  conditions  at  the  first 
slit  were  obtained  from  a-'separate  one-dimensional  calculation  of  a 
charge  in  a rigid  sphere'.  Figure  32  shows  the  pressure  vs.  time  in 
each  chamber  of  the  shield.  Figure  33  shows  the  pressure  distribu- 
tion within  the  shield  0.524  ms  after  shock  impact.  Figure  34  shows 
the  calculated  peak  pressure  just  outside  the  shield.  The  agreement 
with  experimental  data  is  encouraging;  the  peak  pressure  near  the 
shield  would  be  expected  to  be  somewhat  depressed  because  the  pulses 
from  the  exit  slit  have  to  travel  some  distance  before  they  coalesce 
into  a single  sharp  shock  front. 
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5 . 5 Shock  Attenuation  By  The  1/4-Scale  Category  1 I-Beam  Shield 

Figure  35  shows  the  geometry  used  in  the  calculation.  The 
shield  is  a cylinder  with  vertical  rows  of  interleaved  I-beams  forming 
the  outer  wall.  The  slit  widths  along  the  gas  flow  path  are  0.0164, 
0.0128,  0.0154,  0.0204  m;  the  corresponding  vent  fractions  are 
0.217,  0.167,  0.195,  and  0.257. 

5.5.1  Calculations  With  Given  Loading  Pulse 


Figure  36  shows  an  experimentally  determined  loading  pulse 
on  the  inside  of  the  shield.  Figure  36  also  shows  the  calculated 
beam  loads  from  the  shock  flow  thru  the  shield.  Figure  37  shows  the 
calculated  pressures  outside  the  shield. 

5.5.2  Closed  Shield 

Figures  38  thru  42  give  the  calculated  shock  pressures  and 
densities  at  selected  points  on  the  inside  of  the  closed  shield. 

These  rigid-wall  data  are  used  as  inputs  for  calculation  of  the  flow 
thru  the  I-beams. 

5.5.3  Calculations  With  Rigid-Wall  Loading  Pulse 

Figure  43  shows  the  gas  plume  from  an  exit  slit  on  the  side 
of  the  shield.  Figure  44  gives  the  motion  of  the  shock  front  and 
plume  tip.  Figure  45  gives  the  pressure  vs.  distance  profile  at 
1.113  ms.  Figure  46  gives  the  pressure  vs.  time  histories  within  the 
cavities  of  the  shield.  Figure  47  gives  the  peak  pressure  vs. 
distance  outside  the  shield;  the  peak  pressure  is  attenuated  by  a 
factor  of  0.23. 


5.6  Remarks  on  the  Calculation  Method 


The  hydrocode  calculations  have  the  following  disadvantages: 

(a)  The  peaks  of  the  shocks  are  rounded  off  (sometimes  as 
much  as  a factor  of  two)  by  the  artificial  viscosity.  The  true 
waveforms  can  be  estimated  from  a pressure-time  curve  in  the  same 
way  that  an  experimental  gage  record  with  slow  response  can  be  extra- 
polated to  find  the  true  peak.  This  rounding  is  not  a strong 
disadvantage  in  blast  shield  calculations;  it  is  usually  impulse,  not 
peak  pressure,  on  the  shield  elements  that  matters  and  this  is  not 
strongly  affected  by  the  rounding  of  the  peak.  The  most  important 
peak  pressure  is  that  of  the  transmitted  shock  and  this  shock  is  least 
affected  by  rounding  because  it  has  the  slowest  rate  of  decay  behind 
the  front  of  all  the  shocks  in  the  problem. 

(b)  The  artificial  viscosity  introduces  spurious  small 
oscillations.  In  Figures  38  thru  42,  the  underlying  waveform  is 
easily  seen  and  there  is  no  problem  in  ascertaining  the  true  load 
pulse.  In  Figure  30,  however,  the  oscillations  are  mixed  in  with 
the  real  shock  reflections.  They  can  be  separated  out,  if  desired, 
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by  using  smaller  zones  in  the  calculation.  This  reduces  the  period 
of  the  spurious  waves  but  leaves  the  real  waves  unchanged. 

(c)  Costs  can  become  high  if  care  is  not  taken.  The  I-beam 
shield  calculation  shown  in  Figures  43  thru  47  used  165K  (OCTAL)  of 
storage  and  ran  for  25  minutes  on  a CDC  6500  at  a cost  of  $243. 

This  is  judged  a reasonable  cost;  however,  if  the  number  of  zones  in 
each  direction  were  doubled,  the  cost  would  rise  by  a factor  of  ten. 
Some  judgement  is  required  in  setting  up  a problem. 

The  advantages  are  the  ability  to  get  time  dependent  loads 
on  the  shield  elements  and  to  be  able  to  generate  the  complete  flow 
field  so  decisions  can  be  made  on  matters  such  as  where  to  best  place 
the  cavities  and  how  large  to  make  the  vent  areas. 

5.7  Conclusions 


The  calculated  shock  attenuations  are  summarized  in  Table  I. 
The  purely  geometric  vent  fraction,  f,  is  calculated  from 


1/f  = l/f1  + l/f2  + l/f3  + **• 

where  fl,  etc. , are  the  vent  fractions  of  the  individual  plates  or 
gaps  in  the  shield.  The  calculations  show  that  the  shock  attenuation 
by  a given  shield  depends  on  the  character  of  the  shock  that  hits  it. 
The  longer  the  duration  of  the  shock  (i.e.,  the  larger  the  charge), 
the  smaller  the  attenuation.  Large  air  chambers  within  the  shield  are 
advantageous.  Using  a two-dimensional  Euler ian  hydrocode  with  venting 
walls  appears  to  be  a good  way  to  generate  the  detailed  data  needed 
for  improved  shield  design. 
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6.  TIME-DEPENDENT  PRESSURE  LOADS  FOR  BURNING  PROPELLANT  AND 
PYROTECHNICS  IN  A CATEGORY  5 SUPPRESSIVE  SHIELD  STRUCTURE* 

6 . 1 Introduction 

6.1.1  Objective 

The  objective  of  this  task  was  to  determine  the  quantities 
of  rapidly  burning  M10  solid  propellant  and  Mg/NaN03  (45/55) 
illuminant  charge  required  to  load  the  category  5 suppressive  shield 
wall  to  the  proof  test  overpressure  level  produced  by  detonating 
1.11  kg  (2.44  lb)  of  C-4. 

6.1.2  Background 

A category  5 suppressive  shield  is  a vented  enclosure 
designed  to  suppress  an  accidental  detonation  of  an  explosive  igniter 
slurry  mix,  and  to  reduce  the  flame  and  fireball  from  deflagration 
reaction  to  a tolerable  level  (reference  19) . Representative 
operations  to  be  conducted  inside  the  category  5 shield  include: 

(1)  Bulk  propellant  processing 

(2)  Bulk  pyrotechnic  processing 

(3)  Light  metal  or  plastic  HE  (high  explosive)  components 

(4)  Limited  number  of  HE  rounds. 

The  specific  hazard  levels  to  be  satisfied  are: 

(1)  Light  blast  pressure,  less  than  0.34  MPa  (50  psi) 

(2)  Light  to  moderate  fragmentation 

(3)  Flame  propagation  potential. 

6.1.3  Approach 

Computer  program  INBLAS  was  used  to  calculate  the  internal 
overpressure  generated  by  burning  the  propellant/pyrotechnic  inside 
the  suppressive  shield  (reference  15) . The  computer  code  was  modified 
for  this  task  to  handle  finite  burn  rates  in  order  to  allow  burning 
and  venting  to  occur  simultaneously.  The  code  was  also  modified  to 
include  the  presence  of  metallic  additives  such  as  magnesium  and 
sodium  in  the  explosive  compounds.  Of  the  three  hazard  requirements 


♦Symbols  used  in  this  chapter  are  defined  when  first  used. 

19.  Roger,  D.  M.  and  McKown,  G.  L. , "Category  5 Suppressive  Shield 
Test  Report,"  Edgewood  Arsenal  Technical  Report  EM-TR-76001 
(Oct  1975) . 
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(blast  pressure,  fragmentation,  and  flame  propagation)  considered  for 
the  category  5 suppressive  shield,  only  blast  pressure  levels  were 
investigated  here. 

The  weight/burn  area  of  the  propellant  and  pyrotechnic 
burning  inside  the  shield  enclosure  was  varied  in  order  to  determine 
the  quantity  of  each  explosive  composition  which  would  produce  a peak 
internal  overpressure  level  equivalent  to  that  produced  by 
detonating  1.11  kg  (2.44  lb)  of  C-4  within  the  chamber.  The  initial 
overpressure  level  produced  by  detonating  the  C-4  and  the  suppressive 
shield  effective  vent  area  were  also  determined  with  INBLAS;  however, 
in  this  calculation,  instantaneous  burning  (simulating  detonation)  was 
followed  by  venting;  burning  and  venting  did  not  occur  simultaneously. 

6.1.4  Summary 

The  initial  overpressure  level  obtained  by  detonating 
1.11  kg  (2.44  lb)  of  C-4  inside  the  suppressive  shield  enclosure  was 
computed  to  be  0.18  MPa  (26  psi)  for  a 27.4  m3  (966  ft3)  chamber 
volume  with  a vent  area  of  1.603  m2  (17.25  ft2) . 

The  quantity  of  M10  propellant  which  burns  without 
detonating  inside  the  suppressive  shield  enclosure  and  which  generates 
a comparable  internal  blast  pressure  to  that  of  detonating  1.11  kg 
(2.44  lb)  of  C-4  is  17.0  kg  (37.5  lb).  This  result  depends  on  the 
assumption  that  the  total  surface  areas  of  all  M10  propellant  grains 
present  in  the  charge  contribute  to  the  burn  area;  that  is,  all  grain 
surfaces  are  ignited  instantaneously  and  then  burn  uniformly. 

The  peak  overpressure  level  calculated  with  the  magnesium- 
sodium  nitrate  illuminant  mixture  burning  inside  the  supressive  shield 
enclosure  was  found  to  be  a function  of  the  burning  surface  area, 
assumed  constant  during  burning,  and  not  a function  of  the  total 
amount  of  illuminant  mixture  present  in  the  chamber;  however,  a 
minimum  quantity  of  pyrotechnic  must  be  available.  For  example,  with 
a burn  area  of  0.372  m2  (576  in2  - 24  in  x 24  in)  a peak  overpressure 
of  0.1741  MPa  (25.25  psi)  is  obtained  after  having  burned  only  3.27  kg 
(7.22  lb)  of  the  illuminant  mixture.  Therefore,  with  a Mg/NaN03  - 
45/55  weight  somewhat  greater  than  3.27  kg  (7.22  lb)  and  burn  area 
larger  than  0.390  m2  (604  in2)  the  C-4  proof  test  limit  of  0.18  MPa 
(26  psi)  is  exceeded.  However,  if  the  burn  area  is  less  than  0.390  m2 
(604  in2) , then  any  quantity  of  illuminant  mixture,  no  matter  how 
much,  can  be  burned  inside  the  shield  chamber  without  exceeding  the 
proof  test  overpressure  limit. 

6 . 2 INBLAS  Computer  Code 

6.2.1  Original  Version  of  INBLAS 

Program  INBLAS  calculates  the  shock  and  blast  loading 
characteristics  of  an  HE  warhead  detonated  inside  a confining 
structure  with/without  venting;  shock  wave  and  confined-explosion 
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gas  pressure*  loads  are  included.  The  explosives  considered  in  the 
code  are  restricted  to  the  solid  type  of  the  C-H-N-0  form  with 
aluminum  as  the  only  metallic  additive. 

The  code  has  five  options  for  calculating  internal  blast 
phenomena;  three  deal  with  shock  loading  and  two  deal  with  confined- 
explosion  gas  pressure.  The  blast  phenomena  of  interest  in  this 
investigation  are  the  two  confined-explosions  gas  pressure  options 
given  below. 

(1)  Combustion  of  an  explosive  compound  in  a closed 
compartment  and  "instantaneous"  generation  of  the  confined-explosion 
gas  pressure. 

(2)  Combustion  of  an  explosive  compound  and  "instantaneous" 
generation  of  the  confined-explosion  gas  pressure  in  an  initial 
compartment  with  venting  to  additional  vented/unvented  compartments. 

A detailed  description  of  computer  code  INBLAS  including  a 
program  user's  guide  with  a number  of  sample  problems  is  presented  in 
reference  15. 


6.2.2  Modified  Version  of  INBLAS 

Program  INBLAS  has  been  modified  to  consider  a finite 
burning  rate  for  the  explosive  composition,  and  the  code  has  been 
expanded  to  include  metallic  additives  in  the  explosive  composition 
other  than  aluminum,  such  as  magnesium  and  sodium  (reference  20) . 

The  following  data  are  required  for  the  modified  version 
of  INBLAS  if  finite-rate  burning  of  the  explosive  composition  is  to 
be  considered. 

(1)  Table  of  burn  area  versus  explosive  weight  burned. 

The  weight  of  the  metal  additives  is  not  included  in  the  explosive 
weight  (other  than  aluminum) . 

(2)  Table  of  product  of  burn  rate  and  explosive  weight 
density  versus  pressure.  The  weight  of  the  metal  additives  (other 
than  aluminum)  is  not  included  in  the  computation  of  the  explosive 
weight  density. 

If  the  explosive  composition  contains  metallic  additives 
other  than  aluminum  (no  more  than  three  can  be  considered) , then  the 
data  listed  below  must  be  supplied. 


*Conf ined-explosion  gas  pressure  denotes  the  peak  value  of  the  long- 
duration  quasi-static  pressure  which  exists  in  a confining  structure 
or  compartment  following  dissipation  of  the  shock  waves  produced  by 
an  explosion  inside  the  compartment  (reference  15) . 

20.  Lorenz,  R.  A.,  Personal  Communication,  Naval  Surface  Weapons 
Center,  White  Oak  Laboratory,  Silver  Spring,  MD  (Apr  1976). 
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(1)  Weight  of  metal 

(2)  Molecular  weight  of  metal  as  it  appears  in  the  oxide 

(3)  Number  of  moles  of  02  required  to  burn  one  mole  of 
metal 

(4)  Heat  of  formation  for  oxide  ("+"  for  exothermic) 

(5)  Table  of  external  energy  versus  temperature  for  metal* 

(6)  Table  of  internal  energy  versus  temperature  for  metal 
oxide. * 

The  computer  code  was  also  modified  to  include  an  option 
for  considering  a constant  time  step  for  the  burning/venting 
calculation. 

6.3  Model  Description 

6.3.1  Category  5 Suppressive  Shield** 

The  nominal  inside  dimensions  of  the  category  5 suppressive 
shield  cubicle  are  3.20  m (10.5  ft)  by  3.51  m (11.5  ft)  with  a height 
of  2.44  m (8.00  ft)  giving  a volume  of  27.4  m3  (966  ft3) . The  shield 
wall  cross-section  shown  in  Figure  48  (taken  from  reference  21) 
indicates  that  the  panel  structure  consists  of  three  16-gage 
perforated  plates,  four  screens  16  x 16  mesh,  and  two  walls  of 
equally-spaced  angle  lengths.  The  panels,  which  have  an  effective 
venting  area,  were  designed  to  provide  a tortuous  exhaust  path  for 
the  escaping  explosion  product  gases.  The  total  panel  area  is 
35.903  m2  (386.45  ft2) . 

Three  methods  are  discussed  below  for  determining  the 
effective  venting  area  for  the  suppressive  shield  cubicle. 

Method  1.  The  effective  venting  area  (a  = effective  venting  area 
ratio)  is  calculated  for  the  combination  of  angle  stock,  screens,  and 
perforated  plates  shown  in  Figure  48  (reference  21) . 


*The  energy  versus  temperature  data  is  tabulated  in  segments.  Each 
segment  is  defined  by  giving  the  upper  temperature  limit  and 
specifying  the  coefficients  for  up  to  a cubic  fit, 

E (T)  = Aq  + AxT  + A2T2  + A3T3,  to  the  data. 

**  All  dimensions  are  taken  from  reference  21. 

21.  Roger,  D.  M.  and  McKown,  G.  L. , "Category  5 Suppressive  Shield," 
Edgewood  Arsenal  Technical  Memorandum  EM-TM-76001,  EA-4155 
(May  1975) . 
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Effective  Venting  Area  = otEFF  x (Total  Panel  Area) 


Method  2.  The  effective  venting  area  is  computed  using  only  the 
minimum  value  for  a among  the  values  for  the  angle  stock,  screens  and 
perforated  plates  (reference  22) . 


“min  “angle 

Effective  Venting  Area;  = otMIN  x (Total  Panel  Area) 

Method  3.  The  effective  venting  area  is  calculated  with  computer  code 
INBLAS  by  matching  experimental  data  for  C-4  venting  decay  times 
given  in  reference  19  (reference  22) . 

The  first  two  methods  provide  initial  estimates  for 

Method  3. 

6.3.2  Explosive  Compositions 
6. 3. 2.1  Propellant  M10 

Table  II  gives  the  M10  thermodynamic  and  ballistic  data 
which  are  required  input  for  the  burning  option  in  program  INBLAS. 

The  explosive  composition  by  weight  is  determined  from  the  empirical 
formula  (Table  II)  to  be: 

C H N 0 K S 

0.262  0.028  0.129  0.575  0.004  0.002 

Potassium  and  sulfur  are  considered  inert  material. 


22.  Proctor,  J.  F.,  Personal  Communication,  Naval  Surface  Weapons 
Center,  White  Oak  Laboratory,  Silver  Spring,  MD  (Feb  1976). 

23.  Military  Explosives,  Dept,  of  Army  Technical  Manual, 
TM-9-1300-214 , Dept,  of  Air  Force  Technical  Order 
TO-llA-1-34  (Nov  1967) , 
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The  geometry  of  an  M10  propellant  grain  is  defined  in 
Table  III  and  shown  in  Figure  49  (reference  24) . 

The  table  for  burn  area  versus  weight  burned  for 
0.4536  kg  (1.000  lb),  3.086E4  propellant  grains,  is  given  in  Table  IV. 
The  initial  burn  area  is  assumed  to  be  proportional  to  the  total 
weight  of  the  propellant;  for  a charge  weight  of  22.68  kg  (50.00  lb) 
the  initial  burn  area  is  equal  to  62.1  m2  (96200  in2)  which  is  the 
combined  surface  area  of  1.543E6  propellant  grains.  From  the  initial 
point  on  during  the  burning  process  each  grain  present  in  the  charge 
contributes  its  total  surface  area  (as  each  grain  surface  area  varies 
during  burning)  to  the  charge  burn  area;  that  is,  all  grain  surfaces 
are  assumed  to  be  burning. 

6. 3. 2. 2 Illuminant  Mixture 


The  illuminant  mixture  is  composed  of  tumble-mixed 
sodium  nitrate  and  magnesium  granules  (55%  NaN03/45%  Mg) . The 
combustion  products  include  MgO  and  Na20  (reference  19) . 

The  explosive  composition  by  weight  is  given  by: 

N 0 

0.226  0.774 

where  the  weights  of  Mg  and  Na  are  included  in  the  mixture  as  metal 
additives.  For  example,  for  22.68  kg  (50.00  lb)  of  illuminant 
mixture,  the  explosive  weight  is  9.1059  kg  (20.075  lb)  with  10.2  kg 
(22.5  lb)  of  Mg  and  3.368  kg  (7.425  lb)  of  Na. 

Table  V gives  the  illuminant  mixture  thermodynamic  and 
ballistic  data  which  are  required  input  for  INBLAS  (burning  option) . 

Internal  energy/temperature  data  which  must  be  supplied 
to  program  INBLAS  for  metal  additives  (other  than  aluminum)  and  their 
corresponding  oxides  can  be  determined  by  integrating  specific  heat 
(cv)  data  over  temperature.  For  these  calculations,  specific  heat 
data  for  sodium,  sodium  monoxide,  magnesium,  and  magnesium  oxide 
were  taken  from  references  25,  26  and  27.  The  resulting  energy/ 


24.  Military  Standard  Propellants,  Solid,  for  Cannons  Requirements  and 
Packing,  MIL-STD-652C  (MU)  (30  Nov  1973). 

25.  Sittig,  M. , Sodium,  Its  Manufacture,  Properties,  and  Uses, 

Reinhold  Pub.  Corp. , New  York  (1956). 

26.  Touloukian,  Yt  S.  and  Buyco,  E.  H. , Thermophysical  Properties  of 
Matter,  Volume  4:  Specific  Heat,  Metallic  Elements  and  Alloys, 
IFI/Plenum,  New  York  (1970) . " 

27.  Touloukian,  Y.  S.,  Thermophysical  Properties  of  High  Temperature 
Solid  Materials,  Volume  4:  Oxides  and  Their  Solution  and  Mixtures 
- Part  I:  Simple  Oxygen  Compounds  and  Their  Mixtures, 

MacMillan  Co. , New  York  (1967) . 


* 
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temperature  profiles  for  each  of  the  metal  additives  and  their  oxides 
were  then  modified  in  the  vicinity  of  the  melting  temperature  in  order 
to  account  for  the  heat  of  fusion.  This  was  accomplished  by  assuming 
that  the  heat  of  fusion  is  added  to  the  internal  energy  linearly  over 
a temperature  range  of  + 60°K  centered  around  the  melting  point. 

The  burn  area  for  the  illuminant  mixture  is  determined 
by  the  geometry  of  the  storage  container.  The  burn  area  is  assumed 
to  be  constant  during  the  burning  process  and  is  set  equal  to  the 
cross-sectional  area  (top  view)  of  the  container. 

The  burning  rate  for  the  illuminant  mixture  given  in 
Table  V was  determined  from  burning  time  data  (.504  s)  and  burn 
geometry  estimates  (0.406  m x 0.406  m x 0.138  m (16.0  in  x 16.0  in  x 
5.42  in)  - L x W x H)  for  a 22.68  kg  (50.00  lb)  illuminant  mixture 
test  described  in  reference  19. 

6. 3. 2. 3 Composition  C-4 

Table  VI  gives  the  C-4  thermodynamic  data  which  are 
required  input  for  program  INBLAS  (detonation  option) . The  explosive 
composition  by  weight  is  determined  from  the  empirical  formula 
(Table  VI)  to  be: 

C H N O 

0.219  0.036  0.344  0.401 

6 . 4 Results 

6.4.1  Chamber  Venting  Area 

The  venting  area  for  the  suppressive  shield  enclosure  was 
determined  using  program  INBLAS.  Two  estimates  of  the  venting  area 
were  obtained  using  Methods  1.  and  2.  described  in  Section  6.3.1. 

The  results  are  given  below. 

Method  1. 


Effective  Venting  Area  =[  + + J x (Total  Panel 

ANGLE  aPERF  aSCREEN/  Area) 

= 1.005  m2  (10.82  ft2) 


Method  2. 

Effective  Venting  Area  = (aMIN)  x (Total  Panel  Area) 

= 3.590  m2  (38.65  ft2) 
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The  result  determined  by  using  program  INBLAS,  Method  3.,  lies 
between  the  above  two  values. 

Method  3 . 

Effective  Venting  Area  - 1.603  m^  (17.25  ft^) 

Table  VII  gives  a comparison  between  the  C-4  experimental 
data  for  peak  overpressure  and  vent  time  (reference  19)  and  the  INBLAS 
calculations  (using  the  vent  area  value  of  1.603  m2  (17.25  ft2) ) for 
these  variables. 

The  computed  peak  overpressure  value,  PmAX  = 0.18  MPa 
(26  psi)  given  in  Table  VII  for  detonation  of  1.11  kg  (2.44  lb)  of 
C-4  in  the  27.4  m3  (966  ft3)  category  5 chamber  with  a vent  area  of 
1.603  m2  (17.25  ft2) , is  the  proof  test  overpressure  level  for  the 
structure. 

6.4.2  M10  Propellant  Weight 

The  weight  of  rapidly  burning  M10  solid  propellant  which 
would  load  the  category  5 suppressive  shield  wall  to  the  proof  test 
overpressure  level  (vr  0.18  MPa  (26  psi))  was  determined  to  be  17.0  kg 
(37.5  lb) . This  result  is  presented  in  Figure  50  where  peak  over- 
pressure is  given  as  a function  of  propellant  weight.  The  weight  of 
propellant  corresponding  to  a value  of  peak  overpressure  (generated 
during  burning)  in  Figure  50  is  a conservative  estimate;  the  burn 
area  assumption  considers  the  maximum  value  - the  burn  area  is  set 
equal  to  the  total  surface  area  of  all  grains  present  in  the  charge, 
and  all  grain  surface  areas  are  ignited  instantaneously.  In  addition, 
the  burn  is  assumed  to  proceed  in  a uniform  manner. 

The  overpressure  time  histories  for  the  computed  points 
marked  in  Figure  50  are  presented  in  Figures  51(a)  - 51(e).  The 
following  observations  are  made  about  these  computed  profiles. 

(1)  The  peak  overpressure  is  reached  early  in  the  burning  process 
after  only  15-35%  of  the  propellant  has  burned. 

(2)  The  weight  fraction  of  propellant  burned  upon  attaining  the 
overpressure  peak  increases  with  initial  weight  of  propellant. 

(3)  The  peak  overpressure  increases  with  increasing  initial 
weight  of  propellant. 

(4)  The  completion  of  burning  occurs  earlier  in  time  for  larger 
initial  weight  of  propellant. 

(5)  The  free  oxygen  present  in  the  chamber  is  depleted  before 
the  burn  is  completed. 

(6)  The  depletion  of  the  free  oxygen  occurs  earlier  in  time  for 
larger  initial  weight  of  propellant. 
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To  aid  the  discussion  of  the  shape  of  the  overpressure 
profiles  shown  in  Figures  51(a)  - 51(e),  some  details  of  the  computed 
gas  conditions  in  the  category  5 chamber  during  the  4.536  kg 
(10.00  lb)  M10  burn  are  given  in  Table  VIII.  Using  this  table  and 
the  equation  of  state  for  a perfect  gas 

PV  = NRT 

where : 

P = Absolute  pressure 
V = Volume 

N = Number  of  moles  of  gas 
R = Universal  gas  constant 
T = Temperature 

the  peaks  in  each  of  the  overpressure  times  histories  can  be 
explained . 


Since  the  chamber  venting/burning  process  is  accomplished 
with  constant  volume,  the  gas  pressure  depends  on  two  variables,  N 
and  T;  this  dependence  can  be  represented  in  differential  form  by 

dP  _ dN  dT 

P N T 


Table  VIII  lists  the  terms  in  the  above  equation  for  selected  times 
in  venting/burning  sequence  plotted  in  Figure  51(a) 


Referring  to  Table  VIII  and  Figure  51(a),  the  term  — p is 

positive  for  the  20ms  and  40ms  times  (prior  to  peak  overpressure) 

dT 

because  the  gas  temperature  term  — ~ increases  more  rapidly  (from 

1 dN 

energy  released  by  burning)  than  the  gas  density  term  — ^ decreases 

(from  the  net  production/reduction  produced  by  gas  generation/ 
venting)  - hence  the  chamber  pressure  increases.  Beyond  the 

dN 

overpressure  peak,  at  100ms,  the  density  term  — ^ dominates  the 

dT  N 

temperature  term  — in  this  instance,  the  gas  density  reduction 


lowers  the  chamber  pressure  more  severely  than  the  gas  temperature 
increase  can  maintain  or  raise  the  chamber  pressure  - hence,  the 
chamber  pressure  decreases.  Further  along  into  the  burn  sequence,  at 
270ms,  the  chamber  pressure  becomes  relatively  constant  for  a short 
period  of  time  before  the  free  oxygen  supply  is  depleted;  at  the 

dT  dN 

270ms  time,  the  two  terms  — = and  — ^ have  opposite  signs  and  are  almost 
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equal  in  magnitude:*  The  chamber  gas  has  reached  a "pseudo  steady- 
state;"  the  pressure  (energy/volume)  remains  essentially  constant 
because  of  a balance  between  the  energy  being  added  to  the  chamber 
from  burning  and  energy  leaving  the  chamber  by  venting  — during  this 
time,  the  gas  density  continues  to  decrease. 

In  the  M10  burning  calculation,  the  free  oxygen  supply 
within  the  chamber  is  depleted  before  burning  has  been  completed. 

When  the  depletion  of  free  oxygen  occurs  in  the  vicinity  of  the 
"pseudo  steady-state"  flow  regime,  the  chamber  pressure  decays  to  a 
new  "steady-state"  level  after  the  burning/venting  process  has 
achieved  a new  balance  (Figure  51(a))**. 

The  M10  calculation  terminates  following  completion  of 

burning . 

6.4.3  Illuminant  Mixture  Weight 

A quantity  of  illuminant  mixture  in  excess  of  3.27  kg 
(7.22  lb)  with  a burn  area  larger  than  0.390  m2  (604  in2)  will  load 
the  category  5 suppressive  shield  beyond  the  proof  test  overpressure 
level  (v'0.18  MPa  (26  psi)  ) . If  the  burn  area  is  smaller  than  0.390  m2 
(604  in2) , then  any  quantity  of  illuminant  mixture,  no  matter  how 
much,  can  be  burned  inside  the  suppressive  shield  chamber  without 
exceeding  the  proof  test  overpressure  limit.  This  result  is  shown  in 
Figure  52  where  peak  overpressure  is  plotted  versus  burn  area.  The 
peak  overpressure  generated  during  burning  is  a function  of  both 
propellant  weight  and  burn  area;  however,  beyond  a minimum  quantity 
of  illuminant  mixture  (to  provide  sufficient  pyrotechnic  available  to 
allow  the  overpressure  peak  to  be  attained  before  the  burn  has 
terminated)  the  burn  area  is  the  critical  parameter.  All  cases  shown 
in  Figures  52  and  53  were  computed  for  an  initial  composition  weight 
of  22.68  kg  (50.00  lb) . 

The  overpressure  time  histories  for  the  computed  points 
marked  in  Figure  52  are  presented  in  Figures  53(a)  - 53(f).  Comments 
on  these  computed  profiles  are  given  below: 


«*  dT  dN  dP 

*The  method  for  calculating  the  terms  — and  — p (described  in 

Table  VIII)  is  not  accurate  enough  to  provide  good  agreement  between 
dP  dP  dT  dN 

— =r  and  the  sum  — = — = + — r;  when  dp  ^ 0.  This  is  indicated  m 
P P - T N 

Table  VIII  for  the  entries  corresponding  to  270ms. 

**The  exact  time  for  depletion  of  free  oxygen  in  Figure  51 (a) -51(e)  is 
bracketed  by  dashed  vertical  lines  which  mark  adjacent  computer 
printout  steps.  The  shape  of  the  overpressure  time  history  is  denoted 
by  dashed  segments  in  Figures  51(a)  and  (b)  where  abrupt  changes  occur 
because  of  free  oxygen  depletion.  Dashed  profile  segments  are  used  in 
Figures  51(b) -51(e)  for  abrupt  changes  produced  by  completion  of  burn. 
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(1)  The  overpressure  peak  is  reached  early  in  the  burning  process 
after  only  5-25%  of  the  propellant  has  burned. 

(2)  The  weight  fraction  of  propellant  burned  upon  attaining  the 
overpressure  peak  increases  with  burn  area.  The  burn  area  remains 
constant  during  burning. 

(3)  The  peak  overpressure  increases  with  increasing  burn  area. 

(4)  The  completion  of  burning  occurs  earlier  in  time  for  larger 
burn  area.  This  result  is  only  implied  in  Figures  53(a)  - 53(d)  since 
the  complete  burn  is  not  plotted. 

The  expalnation  for  the  appearance  of  the  overpressure  peak 
for  M10  (Figure  51(a)  - 51(e))  given  in  Section  6.4.2  applies  for 
the  illuminant  mixture  results  as  well  (Figures  53(a)  - 53(f)). 

Reference  19  presents  experimental  results  for  rapid 
burning  of  22.68  kg  (50.00  lb)  of  illuminant  mixture  within  a category 
5 suppressive  shield.  Estimates  for  bulk  density,  burning  rate,  and 
storage  container  geometry  were  taken  from  this  test  description  and 
used  for  making  the  computation  of  the  overpressure  time  history  given 
in  Figure  53(b).  The  present  calculation  predicts  an  overpressure 
level  above  .001  MPa  (0.2  psi)  for  a 120ms  duration  with  a 0.08343  MPa 
(12.10  psi)  peak  whereas  reference  19  reports  no  measurable  over- 
pressure during  the  burning  test;  also,  no  measurable  overpressures 
were  detected  during  burning  of  4.54  kg  (10.0  lb)  and  13.61  kg 
(30.00  lb)  of  illuminant  mixture.  Reference  19  states  that  the 
PCB101AOZ  transducers  installed  in  the  shield  for  these  tests  should 
detect  overpressure  levels  above  0.001  MPa  (0.2  psi);  no  additional 
information  about  the  gage,  such  as  full-scale  calibration,  is  given, 
however.  The  same  gages  were  installed  for  C-4  detonation  tests 
where  quasi-static  overpressures  as  high  as  0.076  MPa  (11  psi)  were 
recorded  with  other  gages  (reference  19) . No  measurable  response  was 
detected  with  the  PCB101A02  gage  for  this  overpressure  environment 
either  - however,  this  is  the  calculated  overpressure  level  predicted 
for  the  22.68  kg  (50.00  lb)  illuminant  mixture  experimental  test. 

The  gage  response  in  this  application/calibration  is  questionable. 

Since  the  peak  overpressure  is  overestimated,  the  present 
calculation  for  predicting  the  peak  overpressure  level  attained  by 
burning  the  illuminant  mixture  inside  the  category  5 suppressive 
shield  is  conservative.  Also,  the  values  for  bulk  density,  burning 
rate,  and  proposed  storage  container  geometry  need  to  be  better 
defined  in  order  to  obtain  more  accurate  overpressure  level 
predictions . 
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1 6.5  Comments 

i The  objective  of  this  task  was  accomplished.  The  proof 

test  overpressure  level  produced  by  detonating  1.11  kg  (2.44  lb)  of 
I C-4  in  a category  5 suppressive  shield  was  computed,  and  the 
| quantities  of  M10  solid  propellant  and  Mg/NaN03  (45/55)  illuminant 

l charge  were  determined  which  would  produce  peak  overpressures  during 

\ rapid  burning  at  the  proof  test  overpressure  level. 

Two  different  models  were  used  for ^determining  the  burn 
| surface  area  for  the  explosive  compositions^ 

Model  (1)  - M10  Propellant  - Each  propellant  grain 
contributes  its  total  surface  area  to  the  burning  area,  and  all 
surfaces  are  ignited  instantaneously.  During  burning,  all  propellant 
grains  are  consumed  at  the  same  rate. 

[ Model  (2)  - Illuminant  Mixture  - The  burn  area  is  equated  to 

the  cross-sectional  area  (top  view)  of  the  storage  container.  During 
burning,  the  bulk  pyrotechnic  is  consumed  in  a cigarette  burn  fashion 
from  the  top  to  the  bottom  of  the  container  with  a constant  burn  area. 

\ 

Model  (1)  gives  a very  conservative  estimate  for  the  amount 
of  M10  propellant  required  to  produce  a peak  overpressure  since  the 
maximum  surface  area  available  is  set  equal  to  the  burn  area;  the  peak 
\ overpressures  are  overestimated  for  a given  quantity  of  M10. 

Figure  50  indicates  the  dependence  of  the  peak  overpressure 
on  M10  propellant  weight;  however,  this  figure  can  also  be  re-labelled 
to  represent  the  dependence  of  the  peak  overpressure  to  average  burn 
area.  The  initial  burn  area  (in  the  model)  is  directly  proportional 
to  the  propellant  weight.  Also,  the  burn  area  varies  only 
approximately  15%  during  the  burning  until  all  the  propellant  is 
i consumed.  The  abscissa  for  Figure  50  ranges  from  0-130  m2  in  units 
for  burn  area  (using  an  average  burn  area  for  50  kg  of  M10) . 

Figure  50,  converted  to  burn  area  units,  can  be  used  to  predict  peak 
l overpressures  for  different  estimates  for  burn  area  - providing,  of 
course,  that  the  burn  is  not  completed  before  the  overpressure  peak 
is  attained. 

j; 

The  results  for  Model  (2) , the  burn  area  configuration  for 
the  illuminant  mixture,  depend  on  the  accuracy  of  the  estimates  taken 
for  bulk  density,  burning  rate,  and  storage  container  geometry. 

Figure  52  can  be  used  to  predict  peak  overpressures  for  different 
i estimates  for  burn  rate  and  explosive  bulk  density  if  the  burn  area 

I for  the  illuminant  mixture  is  known.  Since 


w = prA 


I 
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where : 

ft  = Propellant  weight  burn  rate 

p = Explosive  bulk  density  (^40%  of  the  bulk  density  for 
the  composition) 

r = Burn  rate 

A = Burn  area 

then  if  pr  is  changed  to  p'r',  A is  changed  to  A'  such  that  ft  = prA  = 
p'r'A'.  The  new  estimate  for  peak  overpressure  is  read  from  Figure  52 
for  burn  area  A' . 

More  accurate  predictions  of  peak  overpressure  levels  for 
MlO  and  illuminant  mixture  burn  tests  can  be  obtained  if: 

(1)  MlO  — A more  realistic,  less  conservative  model  for 
the  burn  area  is  defined  for  a specific  experimental  test. 

(2)  Illuminant  mixture  — A better  definition  of  bulk 
density,  burning  rate,  and  proposed  storage  container  geometry  is 
established. 
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TABLE  I 

Summary  of  Calculated  Shock  Attenuations 


Type  of 
Shield 

Incident 
Shock 
Overpr . 
P.^  (MPa) 

Transm. 
Shock 
Overpr . 
Pt(MPa) 

0-4 

3.1  step 

1.6 

3.1  decay 

0.085 

T-5 

3.1  step 

0.468 

3.1  decay 

0.096 

BRL  4-plate 
Shock  Tube 
Experiment 

0.483 

decay 

0.074 

0.106 

Transm. 

Factor 

Vpi 

Equiv. 
1-Plate 
Vent  Frac. 

Purely 
Geometric 
Vent  Frac 

0.52 

0.092 

0.025 

0.027 

0.151 

0.0177 

0.029 

0.031 

0.15 

0.22* 

0.017 

0.028* 

0.0195 

♦Underlined  numbers  are  experimental  results. 
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TABLE  II 

Selected  Thermodynamic  and  Ballistic  Properties  for  M10  Propellant3 

Empirical  Formula  (kg  atoms/100  kg) : 

C2 . 1830  H2.7976N0.9180O3.5919S0.0056K0.0113 

Heat  of  Formation,  AHp  (J/kg)  at  298°K:  -2.53 

(cal/g) : -605 

Equivalent  (TNT)  Weight:  0.92 

Propellant  Burning  Rate  m/s:  r = 0.00323  (P (MPa) ) ^ * 695  at  294°K 

from  P=  0 - 75  PMa 

(in/sec):  r = 0.004  (P  (psi)  ) ^ at  70°F 

from  p = 0 - 11,000  psi 

3 

Density  (kg/m  ) : 1670 

(16/in3) : 0.0603 


3 Data  taken  from  reference  23. 
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TABLE  III 

Dimensions  Of  M10  Single  Perforation  Propellant  Grain 

Outside  Diameter  (O.D.  = 3 x W) : 1.41E-3  m (0.0555  in) 

Inside  Diameter  (I.D.):  4.70E-4  m (0.0185  in) 

Length  (L  = 4.5  x O.D.):  6.35E-3  m (0.250  in) 

Web  (W) : 4.70E-4  m (0.0185  in) 


Dimensions  taken  from  reference  24. 
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TABLE  IV 

Burn  Area  Versus  Weight  Burned  For  0.4536  kg  (1.000  lb)  M10 


Weight  Burned  Burn  Area 


kg 

lb 

2 

m 

in2 

0.0 

0.0 

1.241 

1924 

0.02577 

0.05681 

1.232 

1910 

0.05135 

0.1132 

1.223 

1896 

0.07670 

0.1691 

1.214 

1882 

0.1019 

0.2247 

1.205 

1868 

0.1270 

0.2800 

1.196 

1854 

0.1517 

0.3345 

1.187 

1840 

0.1764 

0.3888 

1.178 

1826 

0.2008 

0.4427 

1.169 

1812 

0.2250 

0.4961 

1.159 

1797 

0.2491 

0.5491 

1.150 

1783 

0.2729 

0.6017 

1.141 

1769 

0.2966 

0.6539 

1.132 

1755 

0.3201 

0.7057 

1.123 

1741 

0.3434 

0.7571 

1.114 

1727 

0.3667 

0.8084 

1.105 

1713 

0.3889 

0.8573 

1.097 

1700 

0.4122 

0.9087 

1.087 

1685 

0.4347 

0.9584 

1.078 

1671 

0.4525 

0.9977 

1.070 

1659 

0.4536 

1.0" 

1.070 

1659 
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TABLE  V 

Selected  Thermodynamic  And  Ballistic  Properties 
For  The  Illuminant  Mixture 


NaN03  Heat  of  Formation28,  AHp  (J/kg-N03)  at  298°K:  6.853 

(cal/g-N03) : 1638 

MgO  Heat  of  Formation28 'a,  AHp  (J/kg-mol)  at  298°K:  -601.7 

(kcal/g-mol) : -143.8 

Na20  Heat  of  Formation28 'a,  AHp  (J/kg-mol)  at  298°K:  -504.6 

(kcal/g-mole) : -120.6 


Density  (kg/in  ) : 

1000 

(lb/in8) ; 

0.0361 

Burning  RateC  (m/s) : 

0.28 

(in/s) : 

11. 

28 

Weast,  R.  C. , Handbook  of  Chemistry  and  Physics,  49th  Edition, 
The  Chemical  Rubber  Co.,  Cleveland,  Ohio,  (1968). 

aThe  sign  for  the  heat  of  formation  for  the  metal  oxides  is 
positive  in  computer  code  INBLAS. 

^Estimate  for  tumble-mixed  NaN03~Mg.  The  explosive  density  used 
in  INBLAS  is  40%  of  this  value;  the  weight  of  Mg  and  Na  in  the 
mixture  is  not  included. 

cEstimate  made  using  experimental  data  presented  in  reference  19. 
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TABLE  VI 

3 29 

Selected  Thermodynamic  Properties  For  Composition  C-4  ' 

i 

Empirical  Formula  (kg  Atoms/100  kg)  : C-^  g2H3  54N2  46°2  51 

Heat  of  Formation,  AHF  (J/kg)  at  298°K:  0.139 

(cal/g) : 33.3 

30 

Equivalent  (TNT)  Weight:  1.37  , 


aThe  composition  and  effectiveness  of  composition  C-4  is  unclassified. 
Reference  NAVORD  INSTRUCTION  5511. 4A,  ORD-065,  (29  Dec  1972). 

29 

Dobratz,  B.  M. , "Properties  of  Chemical  Explosives  and  Explosive 
Simulants,"  Lawrence  Livermore  Laboratory  report  UCRL-51319, 

(15  Dec  1972)  . 

30 

Swisdak,  M.  M. , Jr.,  "Explosion  Effects  and  Properties  Part  I - 
Explosion  Effects  in  Air,"  NSWC/WOL/TR  75-116  (6  Oct  1975). 
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TABLE  VII 

Category  5 Suppressive  Shield  Chamber  Venting  Results  For 

Composition  C-4 


Composition  C-4 

Charge  Weight3 
kg (lb) 

PMAX 

Observed3 
MPa  (psi) 

Vent  Time 

Observed3 

ms 

PMAX  h 

Calculated0 
MPa  (psi) 

Vent  Time 

\ 

Calculated 

ms 

0.44  (0.97) 

0.065+0.011 
(9. 4+1. 6) 

40+6 

0.076  (11) 

39 

0.835  (1.84) 

0.123+0.0069 

(17.8+1.0) 

44+2 

0.14  (20) 

48 

1.11  (2.44) 

0.23+0.04 

(33+6) 

38  (?)d 

0.18  (26) 

51 

aData  taken  from  reference  22. 

b 3 3 

Computed  by  program  INBLAS:  Chamber  Volume  = 27.4  m 7(966  ft  )9 

Vent  Area  = 1.603  in  (17.25  ft  ) 

cThe  intercept  for  defining  vent  time  was  determined  by  using  an 
exponential  fit  of  the  form  t = A exp (-Bp)  for  the  last  four  computed 
data  points  (reference  31) . 

^Questionable  Measurement  (reference  22) . 

31 

Ethridge,  N.  H. , "A  Procedure  for  Reading  and  Smoothing  Pressure- 
Time  Data  from  H.E.  and  Nuclear  Explosions,"  Ballistic  Research 
Laboratories  Memorandum  Report  BRL  MR  1691  (Sep  1965) . 
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APPENDIX  A 

Modal  Initial  Velocity  Formulation 

The  methods  developed  for  approximating  solutions  to  plastic 
response  problems  have  all  been  of  the  stationary  mode  type.  A major 
flaw  of  such  methods  is  their  inability  to  correctly  model  the  initial 
conditions  of  motion.  To  illustrate  this  point,  the  case  of  an 
impulsively  loaded  rigid-plastic  beam  with  built-in  ends  will  be 
considered.  The  correct  initial  condition  is 


W(x,o) 


I 

m 


(A-l) 


and  the  exact  solution  to  the  problem  involves  moving  plastic  hinges, 
i.e.,  a non-stationary  mode  solution.  An  excellent  discussion  of  such 
non-stationary  mode  solutions  can  be  found  in  the  report  by  Symonds 
(reference  5) . If  the  values  calculated  by  Equation  (A-l)  were  to  be 
used  as  W(o) , the  modal  initial  velocity,  the  beam  initial  linear 
momentum  and  kinetic  energy  associated  with  the  stationary  mode 
solution  would  be  much  lower  than  the  values  associated  with  the  true 
initial  velocity  field.  Equation  (A-l).  A question  thus  arises  as  to 
what  value  of  W (o)  should  be  used  in  the  stationary  mode  solution. 
Martin  (reference  1)  has  developed  a very  powerful  procedure  for 
choosing  the  "best"  initial  velocity  in  stationary  mode  solutions. 
However,  Martin's  method  applies  only  to  linear  problems  in  which  the 
generalized  stresses  are  constant.  Perrone  (reference  32)  has  pointed 
out  that  in  general  Martin's  method  implies  that  momentum  should  be 
conserved.  If  momentum  is  to  be  conserved,  it  is  found  that 


W (o) 


(A- 2) 


and  for  conservation  of  kinetic  energy. 


- V7  5 


(A-3) 


Numerical  experiments  have  been  performed  using  both  formulations  and 
the  results  indicate  that  Equation  (A-2)  gives  results  in  better 


32.  Perrone,  N. , "Impulsively  Loaded  Strain-Rate-Sensitive  Plates," 
Journal  of  Applied  Mechanics,  Vol.  34,  pp.  380-384  (1967). 
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agreement  with  experiments  than  does  Equation  (A-3) . For  plates,  a 
similar  procedure  yields  the  following  results,  momentum  conserved: 


W (o) 


energy  conserved: 


BL  I 

m[BL/2  - B2  tan  4>/6] 


W(o) 


I BL 

m\[BL/3  - B2  tan  cj>/6] 


(A-4) 


(A-5) 


Again,  as  was  the  case  for  the  beam  analysis,  numerical  experiments 
indicated  that  the  momentum  formulaton  gave  results  more  consistent 
with  experimental  data  than  did  the  energy  formulation. 
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APPENDIX  B 


Numerical  Solution  Procedure 


Of  the  structural  members  considered,  all  except  the  I-beam  were 
shown  to  require  the  solution  of  the  same  basic  set  of  differential 
equations.  However,  as  will  now  be  shown,  the  rectangular  cross 
section  beam  and  the  rectangular  plate  equations  can  be  obtained  as  a 
limiting  case  of  the  I-beam  system  of  differential  equations.  The 
system  of  equations  developed  for  the  I-beam,  Equations  (95)  , (96) , 
and  (97)  can  be  written  in  the  following  general  form. 


Axp(t)  - A2  W = A3  W2  + A4,  0 < W < 6-j^ 

Ajpft)  - A2  W = A5  W2  + A6  W + A?,  6l  < W < <$2 


(Bl) 

(B2) 


A-,p(t) 


a2  W 


Ag  w,  W > §2 . 


(B3) 


If  61  = 62 t the  system  of  Equations  (Bl) , (B2) , and  (B3)  reduces  to  a 
form  identical  to  Equations  (81)  and  (82)  for  rectangular  cross  section 
beams  and  rectangular  plates.  Therefore,  a general  solution  procedure 
for  the  system  given  by  Equations  (Bl) , (B2) , and  (B3)  will  be 
developed  for  use  in  solving  for  the  response  of  any  of  the  structural 
members  previously  discussed. 


Equations  (Bl)  and  (B2)  are  nonlinear,  inhomogeneous,  ordinary 
differential  equations  with  constant  coefficients.  Many  methods  are 
available  for  developing  approximate  solutions  to  such  equations; 
however,  in  choosing  a "best"  method,  one  must  consider  that  the 
solution  should  be  one  which  easily  facilitates  the  calculation  of 
derivatives  of  the  solution  and  times  at  which  the  solution  reaches 
certain  prescribed  values.  It  was  felt  that  the  solution  procedure 
which  best  satisfied  these  criteria  was  a power  series  approximation 
to  the  solution.  To  facilitate  the  evaluation  of  such  a solution,  it 
will  be  assumed  that  the  loading  function  p(t)  can  be  expressed  as  a 
sequence  of  m linear  segments,  i.e.. 


P(t) 


A + b 
n — n 


T , t . < T < t 

n n-1  — n — n 


n 1 ,....,  m 


(B-4) 
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where 


(B-5) 


A power  series  solution  can  now  be  generated  for  Equation  (B2)  subject 
to  loading  of  the  form  given  by  Equation  (B4) . A similar  solution  for 
Equation  (Bl)  may  be  obtained  simply  by  setting  Ag  equal  to  zero  in 

the  solution.  Within  the  interval  t . < t < t , W (t  ) can  be 

n-i  — — n n 

expressed  exactly  by  an  infinite  series. 


M(Tn>  * i-o  Ci  Tn1 


(B-6) 


In  general,  the  series  is  truncated  at  some  point  determined  so  that 
the  solution  has  converged  to  the  desired  accuracy.  For  Equation  (B2) , 
using  eight  terms  in  the  series,  it  was  found  that  the  solution  and 
its  first  derivative  were  accurate  to  a relative  error  on  the  order  of 
10"3.  The  values  of  the  first  eight  coefficients  are. 


Co  = 


Co  = 


C„  = 


Cc  = 


C*  = 


Co  = 


Co  = 

"1‘n-i1 

= w_ 

on 

o 

t— * 

il 

w <Vi) 

• 

' won 

A,-C  2 + A. 
5 o 6 

ci  + A7 

- A1  an 

2A2 

2C0Cl  A5  4 

' C1  A6  - 

A1  bn 

6A2 

-A5(ci2  4 

2 CoC2» 

A6  C2 

12A2 

-V2CoC3 

+ 20^2) 

A6C3 

20A2 

-A5(C22  4 

2C0C4  4 

2c1C3)  - Ag  C, 

30A2 

'A5(2CoC5 

+ 2C^C^ 

+ 2C2C3)  - AgCg 

42A, 


(B-7) 

(B-8) 

(B-9) 

(B-10) 

(B-ll) 

(B-12) 

(B-13) 

(B-14) 
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In  addition  to  this  solution,  one  must  be  obtained  for  Equation  (B3) , 
the  membrane  regime  response  equation.  This  equation  is  an 
inhomogeneous,  harmonic  equation.  If  the  loading  is  assumed  to  be 
expressed  in  a manner  similar  to  Equations  (B4)  and  (B5) , then  the 
solution  for  the  n'th  loading  interval  is, 


W(xn)  = Cg  cos  to  Tn  + Cg  Sin  u)  xn 


(B-15) 


+ 


n 


+ 


where 


(B-16) 


Because  of  the  strong  dependence  of  the  yield  stress  on  the 
strain-rate  for  most  mild  steels,  it  is  not  obvious  what  value  of 
0^  should  be  used  in  calculating  and  M^.  For  bending  only  problems, 

Symonds  (reference  33)  has  proposed  the  use  of  an  average  rotation 
rate  based  on  an  effective  hinge  length.  To  insure  design 
conservatism  in  problems  involving  bending  and  stretching  within 
plastic  hinge  regions,  it  was  decided  in  the  present  formulation  to 

For  a beam  these  can  be 


I,  (B-17) 


II  (B-18) 


rectangular  plates.  The  values 
(B18)  can  be  used  togehter  with 
Equation  (9)  to  calculate  the  dynamic  values  of  and  M^. 

At  this  point,  one  must  ask  whether  a solution  defined  by 
either  Equation  (B6)  or  (B15)  provides  all  the  needed  results,  without 
any  further  computations.  In  general,  the  answer  to  this  question  is 
no,  since  what  is  really  desired  is  the  time  tf  at  which  the  structure 
comes  to  rest  and  the  corresponding  final  displacement,  W(tf) . These 


33.  Symonds,  P.  S.,  "Viscoplastic  Behavior  in  Response  of  Structures 
to  Dynamic  Loading,"  in  Behavior  of  Materials  Under  Dynamic 
Loading , edited  by  N.  J.  Huffington,  Jr.,  ASME  (1965). 


use  average  membrane  strain  rates, 
approximated  by 


e = phase 


e = —j,  phase 
JT 


and  similar  results  are  found  for 
calculated  by  Equations  (B17)  and 
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results  are  not  readily  obvious  from  the  form  of  the  solutions  and  a 
method  for  determining  them  must  be  prescribed.  Mathematically,  the 
following  problems  must  be  addressed: 

1.  solve  W(t*)  = 0,  to  find  out  if  t*  lies  in  the  time 
interval  for  which  the  current  W(t)  is  defined. 

2.  solve  W(t*)  - Si  = 0,  to  find  out  if  the  time  at  which  the 
solution  must  change  to  the  next  region  lies  in  the  current 
time  interval. 

Both  of  these  problems  involve  solving  for  the  root  of  an  equation, 
i .e. , 


F{x)  = 0,  when  x = x*  (B-19) 


A highly  efficient  method  for  solving  such  problems  is  Newton's  method, 
(reference  34) . Newton's  method  involves  iterating  for  x*  by 


F (x*k) 

X*k+1  " x*k  " F' (x*k) 


(B-20) 


Newton's  method  will  always  converge  to  the  smallest  root,  x*  in  the 
region  xo  < x < x*  provided  F' (x)  does  not  have  a zero  in  the  interval 
fx0/  x*]  . First  consider  the  problem  of  finding  a root  of  the 
equation 


F ( t ) = W ( t ) - 5i 


that  lies  in  the  interval  ti  _<  t <_  tf.  In  the  iterative  scheme  let 
= ti  and  suppose  the  desired  root  t*  does  indeed  lie  in  the 
interval.  Newton's  method  will  fail  to  converge  to  t*  only  if 


F' (t)  = W(t)  =0,  ti  < t < t* 


(B-21) 


However,  Equation  (B21)  says  that  the  velocity  is  zero  and  hence  t is 
the  time  at  which  the  structure  comes  to  rest,  i.e.,  the  final  answer 
to  the  problem.  Cases  where  the  iteration  for  roots  of 


F ( t ) = W ( t) 


34.  Isaacson,  E.  and  Keller,  H.  B.,  Analysis  of  Numerical  Methods, 
John  Wiley  and  Sons,  Inc.  (1966) . 
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does  not  converge  are  not  as  easy  to  handle.  If  W(t)  = 0,  for 
ti  < t < tf , the  structure  is  beginning  to  slow  down  at  some  point 
in  the  Interval.  It  is  entirely  possible  that  the  structure  may  come 
to  rest  in  the  interval  t < t < tf.  In  order  to  check  for  such  cases, 
it  is  necessary  to  first  calculate  t and  if  it  lies  in  the  interval 
(ti,  tf) , it  is  then  necessary  to  check  for  roots  of 


W (t)  = 0 


in  the  interval  |t,  tfj  rather  than  jti,  tfj  . 

Computer  programs  which  carry  out  all  of  the  preceding 
caculations  have  been  written  and  the  typical  results  calculated  with 
these  programs  are  given  in  Chapter  4. 
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FIG.  1 . STATIC  COLLAPSE  MODE  FOR  A SIMPLY  SUPPORTED  BEAM. 


FIG.  2.  COMBINED  STATE  OF  GENERALIZED  STRESS  FOR  A RIGID  PLASTIC  BEAM. 
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IMPULSE  (MPa  -s/m) 


FIG.  8.  COMPARISON  OF  I-BEAM  AND  RECTANGULAR  BEAM  LOADED  IMPULSIVELY. 
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FIG.  9.  COMPARISON  OF  WIDE  FLANGE  AND  RECTANGULAR  BEAM  LOADED  IMPULSIVELY  . 
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TEST  CONFIGURATION 

BEAM 
LENGTH 
m (in) 

CHARGE 
WEIGHT 
kg  (lb) 

(r) 

BLAST 

DURATION 

(ms) 

1 

■ ' 

RESPON 

TIME 

(ms) 

3 RING  STIFFENERS 

.762 

7.26 

.32 

0.0143 

.63 

(30) 

(16) 

(.564) 

3 RING  STIFFENERS 

.762 

21.8 

.64 

0.4963 

.95 

(30) 

(48) 

(1.954) 

1 CENTER  RING 

1.52 

7.26 

.32 

.125 

3.16 

STIFFENER 

(60) 

(16) 

(4.95) 

1 CENTER  RING 

1.52 

21.8 

.64 

0.186 

3.46 

STIFFENER 

(60) 

(48) 

(7.31) 

NO  STIFFENER 

3.05 

7.26 

.32 

0.5070 

9.47 

(120) 

(16) 

(19.96) 

NO  STIFFENER 

3.05 

21.8 

.64 

0.5890 

8.81 

(120) 

(48) 

(23.19) 

NOTE:  STRAIN  HARDENING  EFFECTS  ARE  NEGLECTED. 


FIG.  1 1.  COMPARISON  OF  RESULTS  FOR  VARIOUS  CONFIGURATIONS  AND 

CHARGE  WEIGHTS  FOR  THE  UNVENTED  1/4-SCALE  CATEGORY  1 SHIELD. 
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PEAK  PRESSURE  MPa  (psi) 


21.7(3150) 


TEST  CONFIGURATION 

BEAM 
LENGTH 
m (in) 

CHARGE 
WEIGHT 
kg  (lb) 

It) 

BLAST 

DURATION 

(ms) 

FINAL 
DEFL 
m (in) 

RESPONSE 

TIME 

(ms) 

('■ 

};. 

3 RING  STIFFENERS 

0.762 

7.26 

.32 

0.014 

.63 

| 

(30) 

(16) 

(.563) 

3 RING  STIFFENERS 

| 0.762 

| 21.8 

! 

.64 

0.4963 

.95 

i 

(30) 

(48) 

(1.953) 

I 

1 CENTER  RING 

1.52 

7.26 

.32 

0.124 

3.11 

r 

STIFFENER 

(60) 

(16) 

(4.88) 

I1 

1 CENTER  RING 

1.52 

21.8 

.64 

0.184 

3.44 

t; 

STIFFENER 

(60) 

(48) 

(7.26) 

V 

NO  STIFFENER 

3.05 

7.26 

.32 

0.479 

9.04 

I1' 

(120) 

(16) 

(18.87) 

f . 

NO  STIFFENER 

3.05 

21.8 

.64 

0.568 

8.57 

5 

(120) 

(48) 

(22.37) 

NOTE:  STRAIN  HARDENING  EFFECTS  ARE  NEGLECTED 


FIG.  12.  COMPARISON  OF  RESULTS  FOR  VARIOUS  CONFIGURATIONS  AND 
CHARGE  WEIGHTS  FOR  A VENTED  1/4-SCALE  CATEGORY  1 SHIELD. 
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FIG.  13  rNET  FORCE  FROM  A 22  kg  CHARGE  ACTING  ON  THE  WALLS  OF  A 1/4-SCALE  CATEGORY  1 SHIELD. 
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RECORD  FOR  THE  WALLS  OF  THE  1/4-SCALE  SHIELD  UNDER  TRANSIENT  LOAD  . 
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TIME  (ms) 


FIG.  15.  VELOCITY  RECORD  FOR  THE  WALLS  OF  THE  1/4-SCALE  SHIELD  UNDER  TRANSIENT  LOAD, 
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TIME  (ms) 

FIG.17.  VELOCITY  RECORD  FOR  THE  WALLS  OF  THE  1/4-SCALE  SHIELD 
UNDER  TRANSIENT  LOAD  WITH  MULTIPLE  SHOCKS. 
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FIG  18.  DISPLACEMENT  RECORD  FOR  THE  WALLS  OF  THE  1/4-SCALE  SHIELD  UNDER  TRANSIENT  LOAD  WITH 
MULTIPLE  SHOCKS. 
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TRANSMITTED:  TUTTI  10.7  psi  NOTE:  145  psi  = 1 MPa 

EXPT.  15.4  psi 
(REF.  18) 


FIG.  23.  CALCULATION  VS  EXPERIMENT  FOR  BRL  SHOCK  TUBE 
TEST  WITH  FOUR  VENTING  PLATES. 
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ABSOLUTE  PRESSURE  (MPa) 
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t 


TIME  (ms) 

FIG.  26.  PRESSURES  IN  0-4  SHIELD  HIT  BY  A 3.1  MPa  DECAYING  SHOCK  FROM  0.22  kg 
PENTOLITE  CHARGE. 
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ABSOLUTE  PRESSURE  (MPa) 
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PEAK  OVERPRESSURE  (MPa) 
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FIG.  28.  CALCULATED  AND  EXPERIMENTAL  PRESSURES  FOR  0-4  SHIELD. 
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FIG.  29.  ZONING  FOR  CALCULATION  OF  TWO-DIMENSIONAL  FLOW  THRU  T-5  SHIELD 
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FIG.  35,  GEOMETRY  FOR  CALCULATION  OF  22  kg  CHARGE  IN  CLOSED 
CATEGORY  1 QUARTER-SCALE  SHIELD. 
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FIG.  36.  BEAM  LOADS  FOR  1/4-SCALE  CYLINDRICAL  BLAST  SHIELD  WITH  22  kg  CHARGE. 
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DOWNSTREAM  DISTANCE  FROM  SHIELD  (m) 


FIG.  37,  CALCULATED  PRESSURES  OUTSIDE  1/4-SCALE  CYLINDRICAL  BLAST  SHIELD  WITH 
22  kg  CHARGE. 
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FIG.  41,  DENSITY  VS  TIME  AT  CORNER  FOR  22  kg  CHARGE  IN  CATEGORY  1 
QUARTER-SCALE  CLOSED  SHIELD. 
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FIG.  43.  PLUME  FROM  1/4-SCALE  CATEGORY  1 SHIELD  WITH  22  kg  CHARGE. 
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FIG.  44.  MOTION  OF  TRANSMITTED  SHOCK  AND  TIP  OF  PLUME  FOR  1/4-SCALE  CATEGORY  1 
SHIELD  WITH  22  kg  CHARGE. 
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FIG.  46.  PRESSURES  FROM  22  kg  CHARGE  IN  1/4-SCALE  CATEGORY  1 
CYLINDRICAL  SHIELD. 
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4.536  kg  (10.001b)  M10 
27.4  m3  (966  ft3)  VOLUME 


FIG.  51(a).  OVERPRESSURE  TIME  HISTORY  FOR  M10. 
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FIG.  52.  PEAK  OVERPRESSURE  VS  ILLUMINANT  MIXTURE  WEIGHT. 
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FIG.  53(a).  OVERPRESSURE  TIME  HISTORY  FOR  ILLUMINANT  MIXTURE 
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Rock  Island,  Illinois  61201 
Attn:  Mr.  Edward  Jordan 

HQ DA  (DAEN-MCC-I/Mr . L.  Foley) 

Washington,  DC  20314 

HQDA  (DAEN-MCE-D/Mr . R.  Wight) 

Washington,  DC  20314 

Director 

DARCOM  Field  Safety  Activity 
Charlestown,  Indiana  47111 
Attn:  DRXOS-TA/Mr . Olson, 

DRXOS-ES/Mr . Myers, 

TNG+ANAL  DIV/R.  B.  Henry 

Commander 

U.S.  Army  Materiel  Development  & Readiness  Command 
5001  Eisenhower  Avenue 
Alexandria,  Virginia  22333 
Attn : DRCCG 

DRCSF/Mr . W.  Queen 
DRCPM-CS/ COL  Morris 
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Base  Modernization  and  Expansion 

U.S.  Army  Materiel  Development  & Readiness  Command 

Dover,  New  Jersey  07801 

Attn:  DRCPM-PBM-EC/Mr . A.  Dybacki 

Commander 

U.S.  Army  Armament  Command 
Rock  Island,  Illinois  61201 
Attn:  DRSAR-EN/Mr . Ambrosini 

DRSAR- SC/Dr . C.  Hudson 
DRSAR-SF/Mr . R.  Young,  Mr.  J.  Varcho 
DRSAR-PPI/Mr . Serlin,  Mr.  T.  Fetter,  Mr.  S.  Porter 
DRSAR- PP/Mr . A.  Madsen,  Mr.  G.  Cowan,  CPT  Burnsteel 
DRSAR-PPI-A/Mr.  L.  Nelson 

Commander 

DARCOM  Ammunition  Center 
Savanna,  Illinois  61074 
Attn:  SARAC-DE/Mr.  J.  Byrd 

SARAC-DEM/Mr . Huddleston,  Mr.  Sumpterer 

Commander 

Frankford  Arsenal 

Bridge  and  Tacony  Streets 

Philadelphia,  Pennsylvania  19137 

Attn:  Mr.  F.  Fidel,  Mr.  E.  Rempler 

Commander 
Picatinny  Arsenal 
Dover,  New  Jersey  07801 

Attn:  SARPA-MT/Mr.  Saffian,  Mr.  Cana van, 

Mr.  Hickerson,  Mr.  Forsten 


Commander 

USA  Test  and  Evaluation  Command 
Aberdeen  Proving  Ground,  Maryland  21005 
Attn : AMSTE-NB 

Commander 

Dugway  Proving  Ground 
Dugway,  Utah  84022 
Attn:  Dr.  Rothenburg 

Mr.  P.  Miller 


Commander 

Cornhusker  Army  Ammunition  Plant 
Grand  Island,  Nebraska  68801 
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Indiana  Army  Ammunition  Plant 
Charleston,  Indianna  47111 

Commander 

Iowa  Army  Ammunition  Plant 
Burlington,  Iowa  52502 

Commander 

Joliet  Army  Ammunition  Plant 
Joliet,  Illinois  60436 

Commander 

Kansas  Army  Ammunition  Plant 
Parsons,  Kansas  67357 

Commander 

Longhorn  Army  Ammunition  Plant 
Marshall,  Texas  75671 

Commander 

Lone  Star  Army  Ammunition  Plant 
Texarkana,  Texas  75502 

Commander 

Louisiana  Army  Ammunition  Plant 
Shreveport,  Louisianna  71102 

Commander 

Milan  Army  Ammunition  Plant 
Milan,  Tennessee  38358 

Commander 

Radford  Army  Ammunition  Plant 
Radford,  Virginia  24141 

Commander 

Sunflower  Army  Ammunition  Plant 
Lawrence,  Kansas  66044 

Commander 

Lake  City  Army  Ammunition  Plant 
Independence,  Montanna  64056 
Attn:  Mr.  John  Jacobi 

Commander 

Ravenna  Army  Ammunition  Plant 
Ravenna,  Ohio  44266 
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Pine  Bluff  Arsenal 

Pine  Bluff,  Arizona  71611 

Director 

U.S.  Army  Materiel  Systems  Analysis  Activity 
Aberdeen  Proving  Ground,  Maryland  21005 

Director 

U.S.  Army  Ballistics  Research  Laboratories 
Aberdeen  Proving  Ground,  Maryland  21005 
Attn:  Mr.  R.  Vitali 

Division  Engineer  ■■ 

U.S.  Army  Engineer  Division,  Huntsville 
P.O.  Box  1600,  West  Station 
Huntsville,  Alabama  35807 
Attn:  HNDED-CS/Mr.  Dembo 

Mr.  W.  Char 

U.S.  Army  Engineer  Division 
Waterways  Experimental  Station 
P.O.  Box  631 

Vicksburg,  Mississippi  39180 
Director 

DARCOM  Intern  Training  Center 
Red  River  Depot 
Texarkana,  Texas  75502 
Attn:  Dr.  G.  Chiang 

Mr.  W.  H.  Jackson 
Deputy  Manager  for  Engineering 
Atomic  Energy  Commission 
P.O.  Box  E 

Oak  Ridge,  Tennessee  37830 

U.S.  Department  of  Transportation 
Office  of  Hazardous  Materials  Operations 
Transpoint  Building 
2100  2d  Street,  S.W. 

Washington,  DC  20590 

Attn:  Mr.  Erskine  Harton,  Jr. 

Mr.  Frank  Neff 
Mound  Laboratory 
Monsanto  Research  Corp. 

Miamisburg,  Ohio  45342  , 
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Dr.  W.  E.  Baker 
Southwest  Research  Institute 

San  Antonio,  Texas  78284  V 

Division  Engineer 

U.S.  Army  Engineer  Division,  Fort  Belvoir 
Fort  Belvoir,  Virginia  22060 

Mr . John  Komos 
Defense  Supply  Agency 
Cameron  Station 
Alexandria,  Virginia  22314 
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and  Installation  Restoration 
Edgewood  Arsenal 

Aberdeen  Proving  Ground,  Maryland  21010 

Commander 
Edgewood  Arsenal 

Aberdeen  Proving  Ground,  Maryland  21010 
Attn:  SAREA-TD-E 

SAREA-FI 
SAREA-LE 
SAREA-SA 
SAREA-TE 
SAREA-BL-B 
SAREA-BL-E 
SAREA-BL-H 
SAREA-BL-R 
SAREA-BL-T 
SAREA-CL-P 
SAREA-DE-S 
SAREA-MT-TS 
SAREA-MT-TR 
SAREA-MT-H 
SAREA-MT-M 
S ARE A- PA- A 
SAREA-PA-P 
SAREA-PA-Q 

SAREA-TS-R  2 

SAREA-TS-L-  3 

SAREA-TS-E 
SAREA-DM 

* 

Commander 

Aberdeen  Proving  Ground 

Aberdeen  Proving  Ground,  Maryland  21005  * 

Attn:  STEAP-AD-R/RAH 
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Aberdeen  Proving  Ground 

Aberdeen  Proving  Ground,  Maryland  21005 
Attn:  STEAP-TL 

Administrator 

Defense  Documentation  Center 
Cameron  Station 
Alexandria,  Virginia  22314 

Commander 

U.Si  Army  Materiel  & Mechanics  Research  Center 
Watertown,  Massachusetts  02172 
Attn:  DRXMRD/Dr.  C.  Lakshmi-Kantham 

Dr.  Dale  Trott 

Battelle  Columbus  Laboratories 
505  King  Avenue 
Columbus,  Ohio  43201 

Chief  of  Naval  Operations 
Washington,  DC  20350 

Attn:  Operations  Evaluation  Group  (OP  03EG) 

Chief  of  Naval  Research 
800  N.  Quincy  Street 
Arlington,  Virginia  22217 
Attn:  Structural  Mechanics  Branch 

Chief  of  Naval  Material 
Navy  Department 
Washington,  DC  20360 
Attn:  MAT  0323,  Mr.  Jaffe 

Assistant  Secretary  of  the  Navy 
Research  and  Development 
Washington,  DC  20350 

Commanding  Officer 

Naval  Weapons  Evaluation  Facility 

Kirtland,  AFB,  New  Mexico  87117 

Commander 

Naval  Weapons  Center 

China  Lake,  California  93555 
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Commander 

Naval  Ship  Engineering  Center 
Prince  Georges  Center 
Hyattsville,  Maryland  20782 

Commander 

Operational  Test  and  Evaluation  Force 
Naval  Base 

Norfolk,  Virginia  23511 

Director  of  Defense  Research  and  Engineering 
Office  of  the  Secretary  of  Defense 
Washington,  DC  20301 
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